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1. INTRODUCTION 


With the use of composite (non-metallic) materials and 
microelectronics becoming more prevalent in the construction of both 
military and commercial aircraft, the control systems have become more 
susceptible to damage or failure from electromagnetic transients. One 
source of such transients is the lightning discharge. In order to 
study the effects of the lightning discharge on the vital components 
of an aircraft, NASA Langley Research Center undertook a Storm Hazards 
Program in which a specially instrumented F106B jet aircraft was flown 
into active thunderstorms with the intention of being struck by 
lightning. The overall purpose of the program was to enhance the 
capability of detecting and avoiding the hazards associated with 
severe storms and improving design capabilities to protect aircraft 
systems from unavoidable hazards. 

One of the specific purposes of the program was to quantify the 
environmental conditions which are conducive to aircraft lightning 
strikes. To this end, the F106 made over 1400 thunderstorm penetra- 
tions at altitudes ranging from 3.1 km up to 12.2 km (temperatures 
from +5°C to -55°C) and obtained at least 690 direct strikes to the 
airframe (Fisher et al , , 1986). Most of the strikes occurred at 
altitudes in excess of 6 km and temperatures colder than -20°C. Thus, 
a good data base was established for the study of lightning inter- 
action and environment at higher (colder) altitudes (temperatures). 
Analysis of these data has shown that the greatest risk of a lightning 
strike to an airplane occurs at altitudes between 11 and 11.6 km 
(-40 to -45°C) where turbulence and precipitation intensities were 
classified as light (Fisher and Plumer, 1984). Analysis of UHF radar 
observations made during flight operations has indicated that these 
high altitude lightning strikes were triggered by the aircraft 
(Mazur et al . , 1984). 


The results from the Storm Hazards Program are in conflict with 
compilations of data concerning inadvertent lightning strikes to 
aircraft which show that most reported strikes occur in the vicinity 
of 4.5 km altitude, which has been very near the freezing level. 

Fisher and Plumer (1984) also state that the nature of the high alti- 
tude lightning strikes encountered (triggered) by the F106 have been 
generally of a low amplitude current nature, consistent with the known 
characteristics of intracloud lightning discharges (although some of 
the time-resol ved characteristics appear to be different than antici- 
pated). The need for a larger data base for the analysis of lightning 
strike characteristics at lower altitudes where these characteristics 
are anticipated to be somewhat different became evident. The primary 
difficulty during the F106 program was the paucity of lightning 
strikes at altitudes below 6 km. Through 1985, only 75 direct strikes 
to the aircraft were recorded below 6 km, 11% of the total, with 41 
of these 75 occurring in the 1985 field season alone. Virtually all 
of these strikes occurred in the altitude range from 4.3 to 6 km 
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(temperatures colder than 0°C) and were apparently triggered by 
the aircraft (a deduction based on the on-board sensors and 
UHF-band radar data). 

Given the need of obtaining a better data base on lower altitude 
lightning strikes, an effort was undertaken to determine possible 
methods of directing the F106 into lower altitude regions of thunder- 
storms where the probability of a lightning strike would be maximized. 
This effort involved the work of two organizations, the Institute of 
Atmospheric Sciences (IAS) at the South Dakota School of Mines and 
Technology and Electro Magnetic Applications, Inc. (EMA) of Denver, 
Colorado. As the project evolved, additional interests arose in the 
triggering capacity of the aircraft and what effect in-cloud environ- 
mental characteristics had on the triggering phenomenon. Also, the 
characteristics of the lightning discharge itself and its interaction 
with the storm environment were investigated. This report summarizes 
the results of this scientific investigation during its funding period 
from 15 March 1984 through 14 June 1987. 


2. OBJECTIVES 


The objectives of this research program involved the use of the 
IAS two-dimensional. Storm Electrification Model (SEM) to simulate the 
thunderstorm environment in which the F106 was operating. The results 
of the model simulations were initially used in two ways. The first 
use of the model results was to provide EMA with data on the time 
evolution and spatial patterns of the electric fields, small ion con- 
centrations, electric charges on cloud and precipitation particles, 
and the total charge patterns in the simulated cloud. EMA personnel 
then used these data as initial and boundary values in their modeling 
studies of the aircraft lightning environment. 

The main use of the IAS model simulations was to analyze the 
results with an emphasis on looking for relationships between the 
electrical structure of the storm and the basic cloud structure. The 
intent of this analysis was to determine what observable character- 
istics of the cloud correlated with strong electric field regions at 
lower altitudes in the hope of devising a scheme for vectoring the 
F106 into such regions to increase the lightning strike probability. 

An additional objective of the research effort was to develop a 
lightning parameterization scheme for incorporation into the SEM. The 
SEM, in its original configuration, was only capable of simulating the 
development of the electrical aspects of a thunderstorm up to the time 
of first lightning. Without a scheme for simulating lightning, the 
buildup of the electric fields and charges would continue without 
limit, reaching unrealistic values and causing the simulation to 
terminate because of the necessity of reducing the model time step in 
order to maintain numerical stability of the electrical transport 
equations in these high electric fields. In nature, the electrical 
stresses that build up as the charge separation processes proceed, and 
the electric field increases, are relieved by the charge transfer that 
accompanies the lightning discharge. It was felt that an analogue to 
the lightning discharge needed to be incorporated in the model in 
order to accomplish the charge transfer and relief of electric stress 
so that simulations could continue into the mature and decaying stages 
of thunderstorm evolution. Only in this fashion could the complete 
evolution of a thunderstorm be studied dynamically, microphysical ly, 
and electrically. 

In the following sections, the model will be described in detail 
and the results of the various model experiments will be described. 
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3. MODEL DESCRIPTION 


3 . 1 The Base Model 

The theoretical framework is a deep-convection, slab-symmetric, 
two-dimensional, time-dependent (2DTD) cloud model which has been 
applied in the past to several atmospheric convective situations. 
Atmospheric wind, potential temperature, water vapor, cloud liquid and 
ice, rain, snow, and graupel/hail (in the form of ice pellets, frozen 
rain, graupel , and small hail) are the main dependent variables. The 
model has been developed from the works of Orville (1965), Liu and 
Orville (1969), Wisner et al . (1972), Orville and Kopp (1977), and Lin 
et al . (1983). A density-weighted stream function has been used to 
extend the model to deep convection. 

The nonlinear partial differential equations constituting the base 
model include the first and third equations of motion (from which a 
vorticity equation is derived), a thermodynamic equation, and water 
conservation equations (for the three phases). The model has been 
designed such that mesoscale convergence can be superimposed in the 
lower levels and divergence in the upper levels. The manner in which 
such convergence is applied to the model and further details of the 
hydrodynamic equations can be found in Chen and Orville (1980). 

The bulk -water parameterization used in this model divides water 
and ice hydrometeors into five classes: cloud water, cloud ice, rain, 

snow, and high density precipitating ice (graupel/hail) with exponen- 
tial size distributions hypothesized for the three precipitating 
classes. Figure 1 with an accompanying key in Table 1 shows the pri- 
mary cloud microphysical processes simulated in the model. Briefly, 
the production of rain from cloud water is simulated using equations 
based on the works of Kessler (1969) and Berry (1968). Graupel/hail 
is generated by the aggregation of snow, by the capture of snow or 
cloud ice by raindrops (contact freezing), or by the probabilistic 
freezing of raindrops (Bigg, 1953) due to the inherent ice nuclei 
content of a specific volume of water at suitably cold temperatures. 

An approximation to the Bergeron-Findeisen process is used to trans- 
form some of the cloud water to snow. Growth of hail is governed by 
equations for wet and dry growth (Musil, 1970) and shedding of rain 
from hail is included. Cloud water may be transformed to cloud ice 
in the region between 0°C and -40°C using an equation developed by 
Saunders (1%?). Natural cloud ice is normally initiated at tempera- 
tures of -20°C and colder, using an equation after Fletcher (1962) for 
the number of natural ice nuclei active. Homogeneous freezing occurs 
at -40°C. Accretional processes (including riming) involving the 
various forms of liquid and solid hydrometeors are simulated. 

Rain, snow, and graupel/hail have appreciable terminal fall 
velocities, while cloud water and cloud ice have zero terminal velocity 
and, hence, follow the airflow. Evaporation of all forms of hydrometeors 
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TABLE 1 


ORIGINAL PAGE IS 
OF POOR QUALITY 


I 


Symbol 

P IMLT 

P IDW 

P IHOM 

P IACR 

P RACI 

P RAUT 

P RACW 

P REVP 

P RACS 

P SACW 

P SACR 

P SACI 

P SAUT 

P SFW 

P 

SFI 

P SDEP 

P SSUB 

P SMLT 

P GAUT 

P 

r GFR 

P CACW 

P GACI 

P GACR 

P GACS 

P GSUB 

P GMLT 

P GWET 


Key to Figure 1 


Meaning 

Melting of cloud ice to for® cloud wateT, T ^Tq. 

Deposit ion* 1 growth of cloud ice at expense of cloud water. 

Homogeneous freezing of cloud water to for® cloud ice. 

Accretion of rain by cloud ice; produces snow or graupel depending on the amount 
of rain. 

Accretion of cloud ice by rqin; produces snow or graipel depending on the amount 
of rain. 

Autoconversion of cloud water to for® rain. 


Accretion of cloud water by rain. 

Evaporation of rain. 

Accretion of snow by rain; produces graupel if rain or snow exceeds threshold 
and T < Tq. 

Accretion of cloud water by snow; produces snow if T < Tq or rain if T Tq. 

Also enhances snow melting for T T 0 , 

Accretion of rain by snow. For T < Tq, produces graupel if rain or snow exceeds 
threshold; if not, produces snow. For T ^ Tq, the accreted water enhances snow 
melting. 

Accretion of cloud ice by snow. 

Autoconversion (aggregation) of cloud ice to for® snow. 

Bergeron process (deposition and rising) - transfer of cloud water to for® snow. 
Transfer rate of cloud ice to snow through growth of Bergeron process embryos. 
Depositions! growth of snow. 

Sublimation of snow. 


Melting of snow to for® rain, T ^T Q . 

Autoconversion (aggregation) of snow to for® graupel. 
Probabilistic freezing of rain to for® graupel. 
Accretion of cloud water by graupel. 


Accretion of cloud ice by graupel. 


Accretion of rain by graupel. 


Accretion of snow by graupel. 


Sublimation of gravel. 

Melting of graupel to for® rain, T >_ T Q . (In this regime, Pg^ is assuaed to 
be shed as rain.) 


Wet growth of graupel; may involve P^rc and p gaCI miSt * nclu “ c P GACW 01 
p , or both. The amount of wluch is not able to freeze is snea to 

rain. 


or 


i 
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REVP 



Cloud physics processes simulated in the model. See 
Table 1 for an explanation of the symbols. 
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and the melting of snow and graupel/hail are also simulated. For a 
more detailed discussion of this base model, the reader is referred 
particularly to Orville and Kopp (1977) and Lin et al . (1983). 

3.2 The Electrical Model 


This base model has been modified to include various processes 
involved in thundercloud electrification. These modifications 
involved the coupling of the electro-microphysical processes which 
account for the charge exchange between interacting particles within 
the cloud domain. The means of coupling the electro-microphysics to 
the dynamics and microphysics of the base model are based on ideas 
from Chiu (1978), Helsdon (1980), and Kuettner et al . (1981). Each 
of the five classes of hydrometeors is allowed to have a charge asso- 
ciated with it. This charge is then transported with the species of 
hydrometeor according to the general equation 


If ’ - 5 -" ) + k fe (UtQpa) + V-Km7Q * (It) i"ter - < X > 


where Q represents the charge density (C m“ 3 ) carried by the 
hydrometeor class (positive or negative), ^ is the 2D velocity vector, 
Pa is the air density, is the mass-weighted mean terminal fall 
speed of the hydrometeor class, K m is the nonlinear eddy coefficient 
and (sQ/st)i n ter represents the charge exchanged between classes of 
hydrometeors due to interactions. The first term on the right is the 
advection term, the second is the fallout term (zero for cloud water 
and cloud ice charge, which are assumed to move with the airflow), 
and the third is the eddy mixing term. The last term in (1) will be 
discussed below. 

The model framework accounts for the presence of small ions as 
well. The equation governing the number concentration of small ions 
i s 


15.1,2 = -y • [n j 2 V ± n x 2 p i 2 ^ * Km vn i 2 ^ + G - an^ + Src - Sink. (2) 

n L > 5 > 9 


Here n t 2 is the concentration of small ions (the subscript 1 
representing positive ions and 2 the negative ions); y is the small 
ion mobility (m 2 V -1 sec* 1 , pressure dependent); E is the 20 electric 
field vector; G is the ion generation rate due to cosmic ray ioniza- 
tion (ions in* 3 sec* 1 , dependent on height); o is the ion-ion recom- 
bination coefficient (1.6 x 10* 12 m 3 sec* 1 ); Src is the source term 
for small ions from processes other than cosmic ray ionization; and 
Sink is the sink term for small ions from processes other than recom- 
bination. An example of a source process for small ions would be the 
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evaporation of a charged hydrometeor, while a sink process would be 
Wilson capture (ion attachment to hydrometeors). 

With the various types of charges accounted for in the model, a 
total charge density can be diagnosed at each grid point from 


p = e^-nz) + IQ 


( 3 ) 


where p is the total charge density (C m" 3 ) and e is the electronic 
charge (1.6 x 10" 19 C). The first term represents the net charge due 
to the difference in small ion concentrations and the last term repre- 
sents the net charge on the five classes of hydrometeors. With the 
net charge determined, we then use Poisson's equation 


V 2 <J> = p/e 


(4) 


to determine the scalar electric potential, <j> (Volts). In (4), e is 
the electrical permittivity of air (8.86 x 10” 12 kg -1 m" 2 sec 2 C 2 ). 
Finally, the 2D electric field vector can be determined from the 
potential via Gauss's law 



(5) 


The last term in (1) represents the charge exchanged between any 
two classes of hydrometeors during an interaction between those two 
classes. Such interactions may be accretional or non-accretional in 
nature. It should be noted that, in such interactions as P-jacr i n 
Fig. 1 (Table 1), where cloud ice is nucleating rain to form snow, 
the charge on the cloud ice and rain are summed and the total is 
introduced into the snow charge resulting in a three bodied type of 
transfer. 

In the collisional process where a finite separation probability 
exists, two types of electrical interactions are possible: inductive 

processes, where the ambient electric field strength and its orien- 
tation with respect to the particles influence the magnitude and sign 
of the charge transfer; and non-inductive processes, where the ambient 
field exerts no influence, but such factors as temperature, liquid 
water content, and impact speed determine the charge transfer. The 
model has been configured to allow for the accounting of both induc- 
tive and non-inductive processes either individually or in concert. 
This is accomplished by starting with the following representation 
of the charge transferred between two interacting particles. 
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where aQ is the charge transferred due to any non-inductive process, 

Yi is a dimensionless variable which depends on the ratio of the radii 
of the two interacting particles, r|_$ is the vector along the line 
joining the centers of the two particles at the time of impact, r s is 
the radius of the smaller particle, qi_ (q s ) is the charge already 
existing on the large (small) particle, and A (=1-B) is a dimensionless 
variable which depends on the radius ratio of the two hydrometeors. 

The first term, then, represents the non-inductive charge transfer, 
while the second term represents the inductive transfer. The last two 
terms account for the effect of charge already residing on the hydro- 
meteors. Switches are used to activate/deactivate either process. 

A single large hydrometeor may interact with many smaller 
particles in a unit time, so we integrate over the volume swept out 
by the large particle 


9qL = 
at 


/' 


Ef |V S L|N s AqS(a)dA 


(7) 


to obtain the charge acquired per unit time by the large particle. In 
(7), Ef is the collision efficiency between the two particles, |V$l| 
is the relative impact speed, N s is the number concentration of the 
smaller particles, Aq comes from (6), and S(a) is the angle-dependent 
separation probability. The integration results in three possible 
charging rates for the larger hydrometeor depending on which charge 
transfer mechanisms are active 


C i |c 2 + Bq s -Aq|_^j , (8) 

where = Ef |Vsl l^snrL<S> (r|_ being the radius of the large 
particle and <S> is the mean separation probability) and 


-aQ 


C 2 = < 4neyi |E |cos(E ,Vsi_)r s <COSa> 


|E |cos(E,VsL)r s <cosa>-AQ 


non-i nductive, 
inductive. 


( 9 ) 


combined 



Note that in (9) there are two important angles influencing the charge 
exchange: 1) the angle between the electric field vector and relative 

velocity vector of the two particles [cos(E,V$l)]; and 2) the average 
collision contact angle (<cosa>). For all practical purposes, the 
quantity |E | cos(E,Vsl) is adequately approximated by -E z , the 
negative of the vertical electric field component. 

Since (8) implies that multiple charging interactions are possible 
for a precipitating hydrometeor per unit time of fall, and since the 
charge on the large particle appears on the right-hand side of the 
equation, it is not safe, a priori, to assume that 6q|_/<5t is constant 
over a given time step. To account for this, we integrate (8) over 
the length of a time step, At, to obtain the charge transferred to the 
larger particle assuming that all quantities on the right-hand side 
are constant except q|_. This may seem like a contradiction at first 
since both the electric field, |E|, and the charge on the small par- 
ticles, q s , appear on the right-hand side and would seem to be subject 
to time variation as well. It is safe to assume the constancy of q s , 
since it is highly unlikely that any one cloud droplet or ice crystal 
will undergo more than one non-accretional collision per time step. 

In the case of the electric field strength, it is a much more macro- 
scopic quantity, being determined by the entire charged volume of the 
cloud and not so significantly by the local charge changes exhibited 
by individual particles, although the effects of local volumes domi- 
nate. In spite of this, we do admit that the field can change rapidly 
on a local scale, but one other factor exists to help substantiate our 
locally steady assumption. This is the fact that the time step in the 
model is dependent on the electric field strength through the field's 
interaction with the small ion flux term. In order to maintain the 
stability of the ion transport equations, the time step must be 
reduced as the field strength increases. By the time the electric 
field has reached significant levels, the time step has been reduced 
to values of one second or less. Thus, we apply the steady field 
assumption with some confidence. 

Carrying out the integration of (8) results in the following 
expression for the charge accumulated by a precipitating hydrometeor 
during a time step At, 


AqL * (Qm-qLo) , (10) 

where q|_ 0 is the initial charge on the large particle, and the 
parameters Q m and ti are given by 
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Qm 


and 

T i 


(-AQ + Bq s )/A 

> 2 

l (4neYiC0s(E,VLs)r s <C0Sa> + Bq s )/A 

2 

J4neY 1 cos(E,VsL)r's < cosa> - AQ + Bq s )/A 
(E f |V S L|N s nr L 2 <S>A)-i . 


non-inductive, 
inductive, (11) 

combined, 

( 12 ) 


Finally, the charge accumulated by one precipitating hydrometeor 
is multiplied by N|_, the number concentration of the large particles, 
in order to arrive at the total charge density created by interac- 
tions, which appears as part of the last term in (1). The other’ 
components of the last term in (1) include the accretional terms, 
evaporational processes, and ion attachment processes. Examples of 
the formulation of the attachment and accretional processes are 
available in Chiu (1978) and Helsdon (1980). 

3.2.1 Charging processes 

A great deal of laboratory and theoretical work has been, and 
continues to be done to investigate possible mechanisms of charge 
separation. A non-exhausti ve list of possible mechanisms includes: 

1) selective ion capture (Wilson, 1929); 2) convective charging 
(Vonnegut, 1985); 3) induction charging involving non-accretional 
collisions between various classes of particles in the electric field 
(Sartor, 1961, 1967, 1981; Mul ler-Hi 1 lebrand, 1954; Paluch and Sartor, 
1973); 4) effects of freezing potentials during wet growth of hail or 
splashiny collisions between graupel and raindrops (Workman and 
Reynolds, 1948, 1950; Latham and Warwicker, 1980; Shewchuk and 
Iribarne, 1971); 5) non-inductive charging during the riming of 
graupel (Reynolds et al ♦ , 1957; Buser and Aufdermaur, 1977; Takahashi, 
1978; Gaskell and Illingworth, 1980; Ja.yaratne et al. , 1983); 

6) riming of graupel followed by splintering (Latham and Mason, 1961; 
Hallett and Saunders, 1979); and 7) evaporative charging of ice 
crystals in penetrative downdrafts (Telford and Wagner, 1979). It is 
entirely possible that all of these mechanisms (or some as yet 
undiscovered mechanism) are active to one degree or another in various 

thunderstorm situations. However, in order to minimize the degree of 

complexity involved in the modeling work, it was necessary to limit 
the mechanisms included in the model to a workable set. 

The basic charge separation mechanisms included within the model 
can be broken down into four basic categories: 

1) Graupel interacting with particles in a dry growth mode. 

2) Graupel interacting with particles in a wet growth mode. 
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3) Rain interacting with cloud droplets. 

4) Small ions interacting with all classes of hydrometeors. 

Within these various categories, the two types of electrical 
interactions discussed above may be possible, i.e., inductive 
and/or non-inductive transfers. 

The inductive water-water charge transfer interaction is generally 
believed to be ineffective in producing electrification typical of 
thunderstorm conditions on its own. Jennings (1975) has shown that 
such a process is capable of producing early electrification in clouds 
with fields up to 30 kV/m being possible. After this threshold is 
reached, the mechanism extinguishes itself due to the presence of such 
a field causing all collisions to result in permanent coalescence. 

This ability to generate early electrification well below thunderstorm 
values, but significantly above the background fair weather field, may 
affect the nature of other electrification mechanisms which must await 
the formation of ice particles for their initiation. 

The charge separation resulting from the interaction of liquid 
drops and ice particles may be of a non-inductive or inductive nature. 
Kuettner et al . (1981) have summarized the controversy regarding the 
action of a non-inductive, ice-water charging mechanism. They chose 
to include such a process in their model (using conservative values 
for the charge transferred per collision) based on the lack of conclu- 
sive evidence for eliminating the mechanism and the fact that the 
riming process is a well established microphysical phenomenon in 
thunderclouds. 

The non-inductive, ice-ice process is the subject of as much 
controversy as the water-ice mechanism; however, the controversy does 
not center around whether the process is effective, but rather exactly 
what physical mechanism is responsible for the charge exchange and, to 
a lesser degree, the magnitude of the charge transferred. Results 
have been summarized and assessed by Gross (1982). 

The inductive mechanisms involving water-ice and ice-ice 
interactions have been argued as being ineffective in nature because 
only grazing collisions near the electrical equator result in separa- 
tion of particles and charge in the ice-water case, and because of a 
relaxation time limitation on charge migration in the ice-ice case. 

In the ice-water case, it is true that grazing collisions will result 
in minimal charge separation; however. Sartor (1981) has shown that, 
due to the surface roughness of graupel particles, cloud droplets can 
bounce off of such rough particles from almost any impact angle 
negating the above restriction. In addition, the results of Tzur 
and Levin (1981) and Kuettner et al . (1981) have indicated that the 
combination of inductive and non-inductive mechanisms involving the 
ice phase yield the most realistic results. 
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With this discussion in mind and referring to Fig. 1 and Table 1 
for the interactions listed below, the model in its present con- 
figuration can account for both inductive and non-inductive charge 
transfers when rain (Pgacr)» cloud ice (Pgaci)* anc * snow (Pgacs) 
interact with graupel in category (1). In category (2), graupel can 
interact with rain and cloud water (Pgacr an< * Pgacw through Pg We t) 
with the result that excess water is shed as ram (the wet growth 
process). In this category, only inductive charge transfers are 
accounted for due to uncertainty about the actual microphysical 
interactions taking place. In category (3), a limited inductive mode 
transfer is allowed. In addition, both inductive and non-inductive 
charge transfers are allowed when cloud ice interacts with snow. The 
category (4) interactions involving small ions have been described 
above. 

3.3 Initialization 


3.3.1 Environmental initialization 

The base state of the atmosphere for this model is taken from a 
rawinsonde sounding typical of the type of day being studied. This 
sounding of temperature and moisture is, hopefully, representative of 
the atmosphere near the time of the formation of an observed cloud 
which may be compared with the model simulation, if such observations 
are available. The only modification made to the original sounding 
is to make the lowest layer adiabatic (if it is not already in such a 
state). The atmospheric winds are determined by projecting the winds 
from the sounding in the direction of the storm motion, since the model 
is two-dimensional. In addition, the winds determined by this pro- 
cedure are frequently reduced to some percentage of their original 
value because the use of full winds has a tendency to propagate the 
simulated cloud off the grid too rapidly. 

In order to initiate convection, we use a combination of random 
perturbations in temperature and water vapor in the lower 3 km of the 
grid (with maximum amplitude ±0.5°C and ±7%) and a warm bubble in the 
vicinity of the domain center. The bubble has a maximum deviation of 
1.5°C, is 4.8 km wide, and is present between 400 m and 2 km height. 

3.3.2 Electrical initialization 

We assume that an initial steady state exists with ion production 
due to cosmic ray generation, ion loss due to recombination, and ion 
transport due to conduction in the ambient electric field all 
balancing. We ignore ionic diffusion and invoke horizontal homo- 
geneity so that the initial values of the variables are functions of 
height only. Then the aforementioned steady state balance can be 
described by 
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(13a) 


d_ 

dz 


( V* l n 1^2 ) 


G(z) - ctn^ 


and 


(u2 n 2^z) " G ( z ) 


an ^2 


(13b) 


The terms in these equations have been previously defined. E z is 
the assumed fair weather electric field profile patterned after Gish 
(1944) which exhibits an exponential decrease with height such that 


E z = EoCbiexpt-aiZ) + b 2 exp(-a 2 z) + b 3 exp(-a 3 z)] . (14) 


In (14), the parameters E 0 , b-j , and a 3 are varied to adjust the 
vertical profile. This vertical electric field profile can be related 
to the small ion densities by the one-dimensional form of Gauss's law 

dE z= e(n 1 -n 2 ) . (15) 

dz e 

We subtract Eq. (13b) from (13a) and integrate in the vertical to 
obtain 


e ( M 1*1 1 + U2 n 2) " ^T^Z " 3c » 


(16) 


where \j is the total conductivity of the air and J c (A m" 2 ) is the 
fair weather air-earth conduction current, which is assumed constant 
with height. Equations (15) and (16) can be solved simultaneously for 
n x and n 2 yielding 


n i 



1 

e(ui + u 2 ) 


(17a) 


and 
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(17b) 


n 2 = nj 


eE 


z 


e 


If we assume, following Shreve (1970), that the polar ionic mobilities 
vary exponentially with height such that 


mi , 2 = c 1>2 exp(ez) , (18) 

where c t - 1.4 x 10” 4 and c 2 = 1.9 x 10“ 4 m 2 V" 1 sec” 1 , and 
p = 1.4 x 10” 4 m” 1 , then (1/a) and (17b) define the initial profiles 
for positive and negative small ions. Finally, either (13a) or (13b) 
may be solved for the galactic cosmic ray generation function G(z) 
which, after some algebraic manipulation, yields 


G(z) = e(p| + u 2 )^ z ^z + ^z ) 2 + + an x n 2 


(19) 


The initial fair weather electrical state of the atmosphere is 
determined from Eqs. (14) and (17) - (19) by choosing appropriate 
values of E 0 , a-j, b-j , and J c . Since there are generally no specific 
observations upon which to base the selection of these parameters, 
values are chosen that produce profiles which are consistent with 
historical measurements. 

3 . 4 Boundary Conditions 

The top boundary is assumed to be rigid with all variables held 
constant. The vorticity, vertical velocity, rain, snow, graupel, cloud 
water, and cloud ice as well as their associated charges are all set to 
zero. The stream function, entropy, water vapor mixing ratio, small 
ion concentration, and electric potential are maintained undisturbed 
at their initial values. The potential at the top of the model is 
obtained by integrating (5) over the depth of the domain and 
employing (14) for the electric field profile. 

At the lower boundary, the vertical velocity, vorticity, and 
stream function are set to zero. The electric potential is also set 
to zero, consistent with the assumption that the earth's surface is a 
conducting plane. Evaporation and heating rates at the surface are 
prescribed as functions of time. Heat and water vapor are allowed to 
diffuse into the lower boundary. Clouds are not permitted to form at 
the surface, but precipitation is allowed to fall through the surface 
level. Cloud shadow effects (on heating) are also simulated at the 
lower boundary. 
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At the lateral boundaries, the horizontal gradients of the stream 
function and the potential are assumed to be zero. Diffusional 
transport also assumes horizontal gradients of zero for all variables 
at the lateral boundaries. Both inflow and outflow from the domain are 
allowed. The constraints on the potential at the various boundaries 
result in the electric field components E x = 0 at all boundaries and 
E z being calculated at all boundaries. 

3 . 5 Numerical Techniques 

The equations are solved over a 19.2 km x 19.2 km or a 
30 km x 20 km domain with a 200 m grid interval in both the X and Z 
directions. The advection technique used is that of Crowley (1968), 
which is first-order accurate in time and second-order in space. A 
two-step advection scheme is used (Leith, 1965) with vertical advec- 
tion calculated first, followed by horizontal advection. Model 
variables are held constant at lateral inflow boundaries, whereas 
extrapolation via upstream differencing is employed at outflow 
boundaries. The general purpose Helmholtz solver in Cartesian 
coordinates (from the NCAR program library) is used to solve the 
Poisson-type equations for stream function and electrical potential. 
Centered differences are used throughout except at the upper and 
lower boundaries where second order, one-sided differences are 
used. 

3.6 Final Consideration s 

The lightning parameterization has been developed as part of the 
model structure, but will be described later in the context of its 
historical development within the overall framework of this research. 
The following section will enumerate the results of the research. 
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4. RESEARCH RESULTS 


4 . 1 Mai 1 ops Island Simul at ion 

The first use of the SEM in connection with the Storm Hazards 
project was to simulate the 12 September 1983 Operational flight of the 
F106 (flight 83-053) which experienced two lightning events near 28 kft 
during separate cloud penetrations. This simulation is termed a blind 
case because there were no specific observations with which the model 
results could be compared for verification. 

The SEM was initiated using the Wallops Island, Virginia, 00Z 
sounding from 13 September 1983, shown in Fig. 2, which was near the 
time of the flight operations (2215-2258 GMT) and was considered 
characteristic of the atmospheric state in which the storms were 
developing. The temperature has been modified in the lowest hundred 
millibars to make the profile adiabatic, as described in Section 3.3.1. 
Examination of the sounding shows that the Lifted Index is -6.7 and 
the K Index is 42.1, indicating a very high potential for strong 
thunderstorms on that day. In fact, severe thunderstorms might have 
been expected except for the fact that the upper level winds were 
quite light and exhibited little shear. Much convective activity was 
actually observed on that day. With this sounding as input for the 
model, we expected strong convection to develop. 

Figure 3 shows the characteristics of the modeled storm resulting 
from the simulation at two different simulation times: 15 and 22.5 min. 

The left column contains plots at 15 min, while the right column shows 
the same plots at 22.5 min. The plots from top to bottom represent the 
dynamical and microphysical character of the cloud (see plot caption 
for details), radar reflectivity, vertical velocity, total charge 
density, vertical electric field component, and the horizontal electric 
field component. As is evident from Fig. 3, the model simulation of 
this warm-based (~16°C), maritime type storm showed a very rapid 
development with rain formation proceeding via the stochastic coales- 
cence process, well before the ice phase became significant. Note the 
nearly vertical development of the storm due to the lack of wind 
shear. 

The early charge distribution and resulting weak electric fields 
were a result of the limited inductive interactions between raindrops 
and cloud droplets. Once the ice phase appeared, the charge separa- 
tion processes became more vigorous and, in the course of a few 
minutes, charge densities reached the order of 10's of nanocoulombs 
per cubic meter, accompanied by electric field strengths in excess of 
400 kV/m, a value near that necessary to produce lightning breakdown. 

The "classic" thunderstorm dipole structure (positive charge above 
negative) was produced in the simulation with a lower positive charge 
region also in evidence, although there were no observations upon 
which to judge the correctness of these results. 
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Dry adiabats 



Temperature (°C) 

Fig. 2; Wallops Island, VA (WAL) sounding for 00Z, 13 September 1983. 
Solid curve is temperature. Dashed curve is dewpoint. 
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Pressure (mbars) 



H 



Li^l : Characteristics of the Wallops case for two different simulation 

times: left column, 15 min; and right column, 22.5 min. Reading from top 

to bottom, the plots represent: 1) storm dynamical and microphysical 

character including airflow streamlines (dashed lines), cloud boundary 
(solid line), and presence of snow (S), graupel (*), rain (•), and cloud 
ice (-) in amounts greater than a threshold value; 2) radar reflectivity 
in dBz (absolute values systematically too high, but structure is 
representative, contour interval 4 dB at 15 min, 5 dB at 22.5 min; 








Fig. 3 (continued) : 


3) vertical velocity in m sec" 1 (updraft-solid, downdraft-dashed, contour 
interval 3 m sec" 1 ); 4) total charge density in C m - ’ 3 (positive-solid, 
negative-dashed, contour intervals 0.9 pC.m" 3 at 15 min, 1 nC m" 3 at 
22.5 min); 5) vertical electric field in V m" 1 at 15 min, 60 kV m" 1 at 
22.5 min); 6) horizontal electric field (as per the vertical field, 
contour intervals 20 V m" 1 at 15 min; 30 kV m"l at 22.5 min). 









The model run for this Wallops Island case had to be terminated 
after 22.5 minutes of simulated real time due to the continued buildup 
of the electric fields and the resulting prohibitively small time steps 
necessary to maintain computational stability of the ion transport 
equations. A simulation without the electrical processes active was 
run much longer in time. The inclusion of the lightning discharge in 
the model is necessary to overcome this early termination problem. 

Although the electrical simulation had to be terminated 
prematurely, the outputs provided information on the electrical 
environment which might be expected during the early stages of 
thunderstorm development. These data in a form similar to Figs. 3 
and 4, including electric fields, particle densities and charges, and 
ion densities at various altitudes and times within the simulated 
thunderstorm, were sent to EMA for use in conjunction with their 3D 
interaction model (Rudolph et al . , 1984). The plots in Fig. 4 repre- 
sent measurements which an instrumented aircraft would record if it 
made a level pass through the model storm at a given altitude. From 
these plots, EMA scientists were able to extract information con- 
sidered representative of the environment in which the F106 might 
have operated. 

Part of their work involved a determination of the shape factors 
due to the aircraft necessary to interpret the electric field mill data 
which is recorded during each F106 flight. Another aspect of their 
research was to do a parameter study of the aircraft under various 
electrical initial conditions using data from the SEM simulation, where 
appropriate. Their results (Rudolph and Perala, 1985), analyzing the 
data from the electric field mills at the time of a lightning strike to 
the aircraft during flight 83-053, showed that the field was predomi- 
nantly vertical and negative in value. From examination of Fig. 4 at 
22.5 min, the model predicts an environment in the neighborhood of 
28,000 ft which is approaching conditions suitable for lightning 
activity (electric fields in excess of 300 kV/m), indicating that the 
model is capable of generating an environment which is similar to that 
which was observed. 

The investigation of the low altitude strike problem is aided by 
looking at some of the model results. Figure 5 shows contour plots of 
the vertical electric field component (left panel) with solid contours 
indicating an upwardly directed (positive) field and dotted contours 
representing downwardly directed (negative) fields. The interval of 
each contour is 60 kV/m with the zero line indicated. The right panel 
represents the mixing ratio of graupel in the cloud in grams of graupel 
per kilogram of dry air. The contour interval in this panel is 1 g/kg. 
These plots represent the state of these two model variables for the 
Wallops Island simulation at 22.5 minutes after the initiation of the 
run. 


The three horizontal lines drawn through the contour plots 
represent three possible flight altitudes for the F106. The dashed 
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Fig, 4: Examples of measurements which would be recorded by an 

instrumented aircraft flying a horizontal path through the modeled 
storm at 8.4 km AGL and 22.5 min simulation time. Top plot: hail 

(graupel) content in g m~ 3 , number of graupel particles a” 1 with 
diameters > 600 pm, and charge on the graupel particles in nC nr 3 
(positive charge-solid line, negative charge-dashed line). Bottom 
plot: potential above ground in kV, vertical electric field com- 

ponent (E z ) in kV nr 1 (positive-solid, negative-dashed), and 
horizontal field component (E x ) in kV m -1 . 
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In the left panel, downwardly directed fields are indicated by a dotted line and upwardly directed 
fields are contoured with a solid line. The contour interval is 60 kV/m and the zero line is marked. 
In the right panel, the contour lines represent mixing ratios in grams of graupel per kilogram of dry 
air and the contour interval is 1 g/kg. Also represented are possible flight altitudes of the F106 at 






line at 8.4 km represents the altitude at which the aircraft was 
operating at the time of the lightning strikes. The dash-dot line 
represents a flight level of approximately 10 kft, and the solid line, 
a level of around 19 kft. (The major tick marks on the horizontal 
and vertical axes represent distance in kilometers.) These two lower 
altitudes are in the region considered by NASA personnel to constitute 
the low altitude lightning environment (<20 kft). 

Examination of the figure reveals that, at the altitude at which 
the aircraft was operating on that day, a high electric field 
environment would have been encountered (if the model is adequately 
simulating the cloud which was penetrated). The model results further 
indicate that had the F106 been operating in the vicinity of 10 kft, it 
would have encountered relatively low field strengths and probably not 
have been able to trigger a discharge (the triggering process seems to 
predominate at high altitudes). On the other hand, the flight level 
around 19 kft shows that a very concentrated and intense region of 
electric field in excess of 400 kV/m might have been encountered by the 
F106 had it been flying at that altitude. This result seems to indi- 
cate that lower level, high field regions could exist in the types of 
clouds penetrated during the project; however, examination of the 
right panel in the figure shows that this high field region is coin- 
cident with the maximum hai 1 content of the storm (mixing ratio in 
excess of 14 g/kg) . In addition, examination of Fig. 3 shows that 
this high field region also corresponds to a region of strong updraft 
and a sharp updraft/downdraft boundary, indicating the possibility 
of strong turbulence. Therefore, while a volume favorable to the 
initiation of lightning by the F106 in the low altitude region may be 
present, the simulation of this one case indicates that the aircraft 
would have to purposely avoid such a region for safety considerations. 
It is also interesting to note that the majority of the low altitude 
strikes that occurred during the 1985 field season were in the 
15-20 kft region. This leads one to speculate that the low altitude, 
high field region may be a reasonably frequent feature of the storms 
encountered by the F106 and is not always associated with adverse 
flying conditions. More simulations of the electrical evolution of 
such clouds during their mature and dissipating stages would help to 
clarify this question. 

These results brought to light several questions and limitations 
which resulted in additional research. First of all, while the model 
seemed to be producing realistic simulations of the storm conditions 
within which the F106 was operating and the data provided to EMA was 
useful in their analyses, we had no way of knowing whether or not the 
simulations were actually correct. The only way to verify the ability 
of a model to simulate nature is to have actual observations against 
which to compare the model results. The comparison between the F106 
field mill data and the model output done by EMA scientists was 
encouraging, but was not a sufficient test. Such a test has been done 
on the model for a cold-based continental type storm (Helsdon and 
Farley, 198/a, b), giving us confidence that the model can correctly 
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simulate cloud development and its attendant electrification. How- 
ever, there are significant differences between cold-based, conti- 
nental clouds and the warm-based, maritime type clouds found along the 
east coast. These differences are such that we do not feel it is safe 
to assume that, because one type of cloud is well modeled, a different 
type of cloud will be equally well modeled. 

Another, and even more important, consideration involves the 
inability of the model to simulate the electrical development of a 
cloud beyond its early stages. As noted above, the Wallops Island 
simulation had to be terminated after 22.5 minutes simulation time 
because of the constant buildup of the electric field strength. There 
is no mechanism within the model to relax the ever increasing elec- 
trical stresses. Nature takes care of the problem by initiating a 
lightning discharge. In order to be consistent, we realized that it 
was necessary for the model to be able to simulate a lightning dis- 
charge and its effects on the local charges and fields. Only when 
this capability was in place would the model be able to provide simu- 
lations of the mature and dissipating stages of a storm's electrical 
life. 

As a result of the Wallops Island simulation, and through 
discussions with other project personnel, a plan was devised to aid 
in the increase of low altitude strikes to the F106. It was suggested 
that situations be sought where a line of storms was present with 
fairly high reflectivity cores. The F106 should then be directed to 
fly a course parallel to the storm line and upwind of the reflectivity 
maximum in the altitude range from 15 to 20 kft. Based on the Wallops 
simulation, this might put the airplane in the vicinity of the low 
level high field region, but out of the high graupel hazard region. 
Although this recommendation was speculative due to the uncertainty of 
the storm charge structure during the mature stage, it was adopted by 
project personnel . 

4 . 2 Li ghtn i ng Parameteri zati on 

The difficulties with electric field buildup led to investigating 
the question of incorporating (parameterizing) the lightning discharge 
in the SEM. We had two types of lightning to deal with; cloud-to- 
ground and intracloud. While intracloud lightning is the less well 
studied, it is also the most frequent to occur during a storm. It 
also seems, intuitively, to be the less complicated of the two to 
incorporate in the model and therefore was chosen to be examined 
first. 

Without going into the arguments pro and con, we chose to adopt 
the philosophy of Kasemir (1960, 1984) whereby the overall intracloud 
lightning channel is electrically neutral. The charges that exist on 
the channel are the result of ionization along the channel and are not 
"gathered" from the surrounding cloud volume. When one deals with the 
simulation of the lightning discharge, which is a subgrid scale process 
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in the SEM, four basic criteria must be established in order to account 
for the process: 1) initiation; 2) direction of propagation; 

3) termination; and 4) charge redistribution. In attempting a first 
order approximation to the lightning discharge, the idea is to keep 
the processes involved as simple as possible while still trying to 
maintain a physically reasonable approach to the problem. Simplicity 
is also advisable because the physical mechanisms involved in deter- 
mining these criteria are poorly understood at this point. Thus, we 
have basically employed a single parameter approach to the first three 
processes mentioned above. 

Recently Williams et al . (1985) have made a study of electrical 
discharges in polymethlymethacryl ate (PMMA) which had been injected 
with electrons to form various space charge configurations. They found 
adequate scaled correlations to be able to relate the observed behavior 
of the discharges in the PMMA blocks to lightning discharges in thunder 
clouds. Their paper focused on the local electric field strength and 
the space charge distribution as the factors controlling the extent 
and direction of propagation of the discharge channel. While the 
extent and magnitude of the space charge cloud were important in the 
morphology of the laboratory discharges, they identified the local 
electrostatic field as the single most important parameter in con- 
trolling such discharges. And although there is no direct justifi- 
cation for extrapolating their results to what takes place in 
thunderclouds, a scaling approach and the application of their model 
to specific thunderstorm situations argued in favor of accepting a 
possible correlation. As a result of their work and the concept of 
maintaining simplicity in the initial approach to the lightning 
parameterization, it was decided to use the local electric field as 
the parameter controlling the first three criteria in the above list. 

As an initiation criterion, a threshold electric field was chosen 
with a value of 400 kV/m. Williams et al . (1985) quote a value of 
500-1000 kV/m for a breakdown field in their Table 1. It has fre- 
quently been mentioned that the in-cloud breakdown field must be in 
the vicinity of 400 kV/m since the highest reliable in-cloud field 
measurements have peaked near that value. We feel that using this 
value in the model is reasonable because the value at a grid point 
represents the average value of the quantity within a grid box 200 m 
on a side (for this model) and this average would most likely consist 
of both higher and lower values. 

The termination criterion was selected to be a critical field 
strength of 150 kV/m based on work done by Griffiths and Phelps (1976) 
involving the propagation of positive streamers in the laboratory. 

While serious questions can be raised about the appropriateness of 
extrapolating such a laboratory value to the atmosphere, we feel that 
it represents a reasonable threshold to use in the first approximation . 
It is true that the stepped leader which precedes a cloud-to-ground 
discharge appears to propagate through regions where the field strength 
does not exceed a few kilovolts per meter and such an extinction 
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criterion as proposed above would not seem to be applicable; however, 
we are only dealing with the intracloud discharge at present. 

In keeping with our single parameter formulation for the discharge 
process, we have chosen to use the electric field vector in determining 
the direction of the propagation path. Because the model domain is 
composed of grid points in a rectangular pattern, the lightning path is 
constrained to move either along the side of a grid box or along its 
diagonal, depending upon the direction of the electric field vector at 
that point. An algorithm has been developed which computes the angle 
of the field vector with respect to the vertical direction and 
determines whether the path should be taken along the diagonal or 
the side. Even though the vector may be (and usually is) canted with 
respect to the vertical, the diagonal path is not chosen unless the 
canting exceeds a value of 22.5°. Because of these geometric con- 
straints created by the use of a rectangular grid, the resulting 
lightning channel will probably be artificially relegated to a more 
vertical propagation path than would normally be observed in nature; 
however, for a first approximation, this seemed to be a reasonable 
approach which could be improved upon as experience and increased 
knowledge dictated. For example, we recognized that the electric 
field vector at the tip of the propagating leader is important in 
determining the propagation direction; however, this was beyond the 
scope of a first order approximation and was relegated to considera- 
tion as a subsequent improvement. Finally, since the initiation point 
of the channel was the point of highest electric field strength, the 
channel propagation was made bi-directional from that point, moving 
both parallel and anti-parallel to the field vector and terminating 
when the field strength along each segment fell below the cutoff 
threshold. 

In order to evaluate the effectiveness of these criteria, a test 
was conducted using model output from a simulation of the 19 July 1981 
cloud which was observed in the vicinity of Miles City, MT, and pro- 
duced a single intracloud lightning discharge (inferred from the 
electric field measurements of the NCAR/NOAA sailplane; see Dye et 
al . , 1986). The initiation, propagation, and termination criteria as 
outlined above were applied to the model predicted state of the cloud 
at a time when the electric field was sufficient at some point in the 
domain to exceed the initiation threshold. From this point, the 
discharge path was computed in two directions following the electric 
field vector. The upward and downward paths continued until their 
termination thresholds were reached, at which point the discharge in 
that direction was stopped. 

The results of this process are shown in Fig. 6 which depicts the 
model calculated lightning channel superimposed on the ambient charge 
distribution (left panel) and the vertical electric field component 
(right panel). We can see that the discharge originates in the region 
of maximum vertical electric field, which corresponds to a region of 
low net charge density. The discharge path is basically vertical 
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following the dominant vertical electric field component and 
penetrates both the positive (upper) and negative (lower) charge 
centers. One should note, however, that the effect of the horizontal 
component of the electric field vector (not shown) becomes sufficient 
in the vicinity of the downward propagating part of the path so that 
the channel is deflected to the left as termination takes place. The 
total path length of the discharge channel in this instance was 
2682 meters. 

The parameterization of the lightning process was completed by 
specifying the manner in which the charge created along the channel 
is redistributed into the environment. This aspect of the problem was 
approached in the following manner. Using ideas from Kasemir (1960, 
1984), we assume that the lightning channel is going to be elec- 
trically neutral over its entire length, indicating the deposition of 
equal amounts of positive and negative charge along the two portions 
of the simulated discharge. In the model, this charge deposition is 
accomplished by assuming that the charge density per unit length along 
the channel is proportional to the difference between the ambient 
potential, <j> (in Volts), at the grid point in question and the 
potential, <j, 0 , at the initiation point of the discharge. Thus 


Qo = -k (♦*♦<>) 


( 20 ) 


where Q 0 is the linear charge density at a point along the channel, 
and k is a proportionality constant (farads/m). The negative sign 
assures that the charge deposited along that portion of the channel 
is basically opposite in sign to the ambient net charge in that grid 
volume, since the sign of the potential and the charge responsible 
for that potential are generally the same. 

The application of Eq. (20) results in a charge distribution 
along the channel similar to that shown in Fig. 7 between +1.2 and 
-1.5 km. The two tails extending beyond these points in Fig. 7 will 
be explained subsequently. Since we are working with a discrete grid 
and the termination criterion forces the channel path to stop at a 
specific grid point (not at the exact place where this criterion is 
just met), charge neutrality which is required by theory is not 
guaranteed by this process. In order to overcome this difficulty, 
we arbitrarily extend the channel path four grid points beyond the 
designated termination point at each end. We first integrate over the 
channel path to find the total charge deposited on the sections of 
opposite polarity. We then compare these two and, at the end of the 
path segment having the greatest total charge magnitude, we cause an 
exponentially decreasing amount of charge to be deposited over the 
four additional grid points using the following expression 


Qj = Q n exp(-k{[j-n]A} 2 ) 


( 21 ) 
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ng the channel path using Eq. (20) and a value of 
channel is delineated by the vertical bars on 




where j = n±l,...,n±4, Q n is the linear charge density (C/m) at the 
termination point of that segment [as determined by Eq. (20)], a is 
the grid spacing, and k is chosen such that Q n± 4 = Q n /1000. This new 
segment of charge density is integrated to obtain the total charge 
deposited along the segment of that polarity. We then calculate the 
absolute difference between the charge residing on the completed 
segment and the charge residing on the segment of opposite polarity. 
This difference is used as a basis for calculating a charge density 
decrease over the four yrid points at the end of the remaining 
segment. 

As with the charge calculations on the initial channel, a 
trapezoidal integration is performed using Eq. (21) to represent the 
charge decrease over the four grid points. An iterative process 
(Newton-Raphson) is used on the integration scheme with the exponen- 
tial factor k being the unknown quantity. By varying k, the rate of 
decrease of charge density away from the initial endpoint of the 
segment is controlled and the total additional charge distributed off 
the end of the segment can be set so that the total charge of that 
segment just balances the total charge on the other segment. By this 
procedure, the neutrality of the discharge channel is obtained over 
the discrete grid. 

Figure 7 shows the linear charge density calculated along the 
channel using Eq. (20) and a value of k = 2.5 x 10 -11 farads/m. The 
original channel is delineated by the 2 vertical bars on the hori- 
zontal axis. The effect of the charge balancing procedure is evident 
at the 2 tails. The total channel length is nearly 4.5 km. For this 
particular discharge (value of k), the total charge transferred is 
12 C. This may be a bit high for an intracloud discharge; however, 
it is not an unreasonable value. 

To this point, we have the lightning discharge represented as a 
path on the discrete grid with a linear charge density at every grid 
point along the path. In order to allow this discharge to interact 
with the other model parameters, we need to make some further modifi- 
cations. The discrete discharge path represents a line (ID) source in 
the 2D model. Because of the numerical methods employed in the model 
calculations, the sudden appearance of a line source of ions would 
represent a shock which the model could not accommodate. In order to 
make the channel ion production tractable for the model numerics, the 
charge must be distributed to some of the grid points adjacent to the 
channel path. A minimum of three grid points must be involved for 
numerical stability, the channel path itself and one grid point on 
each side. We have chosen to employ a total of 9 grid points (the 
path grid point plus four on each side) to represent the charge 
distribution. An expression identical to Eq. (21) is used to specify 
the decay of the magnitude as a function of distance away from the 
channel, again with Q n± 4 = Q n / 1000 . For either a vertical or diagonal 
path orientation, the charge spreading is done in a horizontal 
di rection. 
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One additional consideration must be made. The charge density 
along the channel is expressed in terms of C/m. This must ultimately 
be converted into an ion density (ions/m 3 ). The linear charge den- 
sity, expressed by Q n , represents the charge density through a grid 
volume which surrounds the grid point and, therefore, can be converted 
to an equivalent volumetric charge density, assuming that the charge 
is distributed uniformly throughout the grid volume in question. In 
addition, when we undertake to distribute the line charge laterally 
away from the discharge path, we would artificially create charge 
i f we do not take steps to maintain charge conservation. However, 
charge conservation can be accomplished during the spreading by 
adopting the following procedure. 

Figure 8 is a representation of the charge conservation concept. 

Q n represents the linear charge density at a particular grid point 
along the discharge path. Assuming that Q n represents the charge 
distribution in the grid volume AXAyAz surrounding the grid point, 
then the total charge in the grid volume is 


QT = QnAl 


( 22 ) 


where a1 is the length of the discharge path segment for that grid 
volume. Because we are using a two-dimensional model, we must 
generate a fictitious volume with which to work. Although, in a 
strict sense, the adoption of slab symmetry for the 2D model implies 
an infinite extension to all model features in the y-direction, we 
choose here to adopt a concept of infinity which encompasses only the 
range of influence of the process being considered. Since we are 
employing a 9 grid point spreading procedure in the x-direction, we 
hypothesize a 9 grid point spreading in the y-direction also to 
comprise the influence volume. 

In Fig. 8, the rectangular volume at the center represents the 
total charge associated with the grid volume in question. The two 
Gaussian type curves represent the redistribution of that charge over 
the 9x9 rectangular grid array surrounding the grid point where the 
total charge in both distributions is the same. 

Mathematically the procedure is developed as follows. The total 
charge under the two distributions is equated such that 


q T - Q n Ai 


Al 4 Ay 4ax _ k 2 . 2 

= Q'J J J e * lX e k2 ^ dxdydz 


o -4Ay -4 ax 


which becomes 


(23) 
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defined by the Gaussian curves are equal. 
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(24) 
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which can be separated to yield 


4 *y _k v 2 4ax _b y 2 
Q t = Q Al = 4Q*Alf e k 2 y dy r e k i x dx 

' J o J o 


(25) 


Here, Q' is the magnitude of the charge density (C/m 3 ) at the grid 
point under consideration which will result in the same total charge, 
Qy, when distributed under the exponential distribution rule. 

The two integrals are evaluated using 6 point Gauss-Legendre 
Quadrature with kj = k 2 , yielding a fall off to 1/1000 of the central 
point value at the 4th grid point [see discussion following Eq. (21)]. 
The result of this integration yields a relationship between Q* and Q n 
such that Q' = Q n /2 . 909 x 10 5 . Therefore, the linear charge density 
at a grid point may easily be converted to a volumetric charge density 
spread out laterally away from the discharge path in a manner identi- 
cal to Eq. (21), but with Q' replacing Q n . Figure 9 shows the results 
of the horizontal spreading of the charge distribution depicted in 
Fig. 8. because of the exponential fall off, most of the charge is 
concentrated nearest to the main channel. The upper portion of the 
channel contains negative charge and the lower portion is positive. 

This process conserves the original total charge transferred by the 
discharge. 

Each calculated grid point charge density must be converted to an 
equivalent ion density. This is accomplished by dividing each charge 
magnitude by the electronic charge, e = 1.6 x 10 -19 C/electron, 
assuming that all the ions produced are singly charged. The ions 
resulting from the redistribution process are then added into the 
ambient ion fields as they exist at that point in the simulation. 

With the ion fields modified, a new total charge density is calcu- 
lated, then Poisson's equation is employed to diagnose a new potential 
field from which the modified electric field components are calculated. 
This completes the discharge process in the model and the time marching 
is resumed. The ions which were injected into the cloud along the 
discharge path interact with the charged hydrometeors in subsequent 
time steps and continue to modify the charge distribution and the 
resulting electric fields. 


! 
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4 . 3 Experiments with the Li ghtni ng Parameterization 


4.3.1 Charge transfer parameter study 

The first experiment performed with the parameterization scheme 
used the value of k referred to above and resulted in a very strong 
influence on the subsequent electrification. Because of this, we 
decided to try a series of k values to see how differing amounts of 
charge transfer affected the electrical behavior. Using 10 values of 
k ranging from 2.5 x 10” 12 (weak charge transfer) to 2.5 x 10" 11 
(strong charge transfer), we found a vastly different character to the 
electrical structure of the cloud. These results ranged from a minor 
perturbation for the weak cases to a complete domination of the charge 
structure and electrical evolution for the strongest cases. The 
19 July CCUPE simulation was used as the framework for these 
experiments . 

For demonstration purposes, some aspects of four of the cases 
have been chosen to illustrate these differences. In the following 
graphs, the numerical labels represent the following values of k [in 
Eq. (20)] and amounts of charge transferred: 

1 -> k = 2.5 x 10-n, 12 Coul 

2 -> k = 2.5 x 10“ 12 , 1.2 Coul 

5 -> k = 1.0 x 10-n, 4.8 coul 

7 -> k = 1.5 x 10-n, 7.2 Coul. 

Figure 10 shows the variation of the maximum snow charge 
density (a), the maximum negative vertical (b), and maximum negative 
horizontal (c) electric field strengths for a period of 30 seconds 
following the time of the simulated lightning discharge for each of 
the four cases listed above. The majority of the snow in the simu- 
lated cloud is charged positively as a result of riming electrifica- 
tion. The lightning discharge occurs in the region of maximum snow 
charge density with both positive and negative ions being deposited 
amongst the positively charged snow. The negative lightning ions are 
formed in the region of maximally charged snow. Figure 11 shows the 
snow charge density before and after the discharge for case 5. 

Figure 10a shows that, in all four cases, the negative ions attach to 
the snow particles to reduce the maximum charge density. For the 
weakest case (2), the effect is very small and the electrical pro- 
cesses tend to return to their previous state. For the intermediate 
cases (5 and 7), there are sufficient negative ions for the attachment 
process to continue to erode the charge on the snow. For the 
strongest case (1), the trend is the same initially, but then shows a 
large increase in the snow charge. Because the data used to make up 
this plot are the domain maximum values, this large increase repre- 
sents the effect of the attachment of positive ions from the lower 
portion of the channel to snow particles lower in the cloud. While 
the original maximum charge region has been eroded by the negative 
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Fig. 11 : Snow charge density at 58.5 min (before lightning) and 

59.0 mi n (30 sec after lightning). Solid lines represent positive 
charge and dashed lines, negative charge. The attachment of lightning 
produces ions (both positive and negative) is apparent. 
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ions, a new region of positive maximum has been created. Figure 11 
shows the results of the discharge ion injection on the snow charge 
after 30 seconds. 

Figure 10b shows the effect of the discharge ions on the maximum 
negative vertical electric field component. The weak and intermediate 
discharges all exhibit a tendency to reduce the field strength, with 
stronger discharges exhibiting the stronger tendency. This seems 
intuitively satisfying because we tend to think of lightning as acting 
to reduce the electrical stress in a cloud. For the strong case, how- 
ever, the results are somewhat unexpected. Here there is an initial 
i ncrease in the magnitude of the negative field component followed by 
a return to a near steady state at almost the initial value. Examina- 
tion of the structure of the vertical field component for cases 2, 5, 
and 7 shows that the lightning causes an ever increasing perturbation 
on that structure; however, for case 1 an entirely different structure 
is established by the discharge. 

Figure 10c, which shows the maximum negative horizontal electric 
field component, further demonstrates the transition from perturbation 
to domination. Here, only the weak case shows a decrease in the field 
strength with a rapid return to the slow field buildup. The two 
intermediate cases both show an increase in the field strength imme- 
diately following the discharge; however, case 5 then shows a field 
reduction below the initial value, while case 7 continues to maintain 
field strengths greater than the initial value. The strong case shows 
a dramatic increase in this field component and maintains the 
increased magnitudes. 

4.3.2 Comparison with lightning induced electric field changes 

These widely varying results led us to conclude that there was a 
great deal that we did not understand about how the lightning discharge 
operates in nature. Therefore, we saw that there was a need to 
continue to study how the parameterization scheme worked within the 
context of the model and what the various results meant relative to an 
actual thunderstorm. Unfortunately, in order to interpret our results 
as to their meaning within the realm of actual storms, we needed to 
have actual observations for comparison. This is one of the problems 
confronting our investigations. There are very few actual in situ 
observations of lightning and its effects on the electrical structure 
of thunderstorms. Most of the measurements that are available are of 
the change in the electric field associated with lightning made by 
suitably equipped balloons or aircraft which happen to be flying in the 
vicinity of discharge. In the case of our use of the 19 July CCOPE 
case as a framework within which to develop the lightning 
parameterization, one such observation was available. 

During the CCOPE experiment, the National Center for Atmospheric 
Research flew a sailplane equipped with a cylindrical field mill for 
measuring two components of the electric field. This sailplane took 
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measurements during an ascent in the updraft of the 19 July cloud and 
during two passes along the periphery of the cloud after exiting the 
updraft. One pass was made along the edge of the cloud at an altitude 
of b.3 km MSI. (4.5 km AGL in the model) and detected an apparent 
intracloud discharge. The electric field recorded by the sailplane 
during its pass along the cloud edge is shown in Fig. 12, which is 
reproduced from Fig. 9 of Dye et al . (1986). The sharp discontinuity 
in the horizontal field component just after 1637 MDT led them to infer 
the occurrence of an intracloud discharge. 

In order to compare this observation with the model simulation, we 
had to establish a point in the model domain to use as a reference for 
extracting the electric field data. We assumed that the sailplane 
would fly perpendicular to the model plane (X-Z plane). We used the 
simulated cloud boundary for the horizontal position and an altitude of 

4.4 km AGL for the vertical position of the sailplane. Just prior to 
the lightning discharge in the model, the sailplane would measure a 
horizontal electric field component of about 6b kV/m. The actual 
horizontal component measured by the sailplane just prior to the 
discharge was 7 kV/m. While the modeled and observed values differ by 
an order of magnitude, the horizontal position of the model sampling 
point was quite arbitrary and the field strengths vary rapidly in the 
horizontal. Thus, the discrepancy in the quantitative aspects of the 
comparison is not particularly significant. After the simulated 
intracloud discharge, the sailplane would measure a horizontal field 
component value of about 55 kV/m. The field change due to the 
simulated intracloud lightning was -10 kV/m. The actual horizontal 
electric field measured by the sailplane following the discharge was 

2.5 kV/m, yielding an actual field change due to lightning of about 
-4.5 kV/m. 

Figure 13 illustrates the time evolution of the modeled horizontal 
and vertical electric field components as might be measured by the 
sailplane at the time of the intracloud discharge. We can compare 
these with the actual curves shown in Fig. 12. Allowing for the 
difference in magnitudes, there is a high degree of similarity between 
the two curves representing the horizontal component. But we also 
note that, for the vertical field component, the model result does not 
compare well with the actual sailplane data. This discrepancy may be 
the result of model geometry. Since, in this two-dimensional model, we 
assume that the cloud is homogeneous and infinite in the Y direction, 
the cloud geometries are not equivalent. On the positive side we note 
that, in thunderclouds studied in New Mexico, balloon borne field mills 
have recorded simultaneous discontinuities similar to those in Fig. 13 
for all three components of the electric field (cf. Weber et al . , 

1982). 

4.3.3 Lightning ion interaction with cloud particles 

In addition to the study of the effect of the lightning discharge 
on the electric field structure, another point of interest was how the 
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Fig. 12 : Reflectivity contours near the time of the intracloud 

discharge with tracks of the NCAR sailplane and the Desert Research 
Institute Aerocommander superimposed. Arrowheads show the position of 
the aircraft on each minute. Sailplane measurements of electric field 
are shown above. The discontinuity in the horizontal field trace is 
evident. (Extracted from Dye et_ aj_. , 1986, with permission). 
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Fig, 13 ; Time evolution of the horizontal (A) and vertical (B) 
electric field components at the assumed sailplane position 
around the time of the lightning discharge. 
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ions produced by the lightning interacted with the hydrometeors in the 
cloud. This is an area in which little in the way of observational 
data are available, but which lends itself nicely to study within the 
model. To do this, we simulated a number of aircraft traverses through 
the storm at an altitude of z=4 km (AGL) and at 7.5 second time 
intervals. The 4 km height was near the lower termination point of the 
discharge channel. The results concerning the charge on the cloud 
droplets are shown in Fig. 14, while the charge on graupel particles 
is shown in Fig. 15. 

Before the lightning, the cloud water charge has two positively 
charged regions (one located at x=13 km and the other at x=17 km) and 
two regions of negative charge centered at x=14 and x=16.5 km, 
respectively. For the graupel, there is a broad region of negative 
charge centered near x=16 km. The lightning discharge produces a large 
number of ions (depending on the amount of charge transferred by the 
channel) which can interact with hydrometeors in the vicinity by the 
ion attachment process (Wilson capture). In Fig. 14, we see that the 
region of strong negative cloud water charge density near the lightning 
channel at x=17 km is reduced rapidly, disappearing after 58.625 min, 
while the adjacent region of positive charge broadens and intensifies 
with time. Referring to Fig. 15, we see that the lightning induced 
change in the graupel charge distribution is initially apparent as a 
notch developing in the negative region at x=17 km, the channel axis. 

We note that the effect is small until the negative cloud water region 
has disappeared. After that time, the positive ions which were 
produced by the lightning channel interact more strongly with the 
negative graupel charge. At 59 min, the graupel charge distribution is 
characterized by two peaks. The adjacent region of positive cloud 
water charge becomes stronger during this time period, but more slowly. 
The cloud droplets are a more efficient sink for ions by the attachment 
process than are the graupel particles. This is because the total 
surface area associated with the cloud droplets is usually much larger 
than that associated with graupel even when the cloud water content is 
smaller than the water content contained in the graupel. The ion 
attachment process is dependent on the surface area of the attaching 
species. 

4.3.4 Energy dissipation 

Another way to check the consistency of the lightning 
parameterization is to look at the energy associated with the 
discharge. It is widely accepted that the energy for the lightning 
discharge is provided by the ambient electric field of the thundercloud 
environment as well as the electric field at the tip of the leader. 

The energy which is dissipated by the lightning is used for ionizing 
and heating the channel, producing electromagnetic radiation, building 
a magnetic field, producing acoustical waves, etc. In order for the 
parameterization scheme developed for the SEM to be valid, the 
electrical energy dissipated by the model channel should agree in order 
of magnitude with observed or theoretically estimated discharge energy 
dissipation. 
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Fig. 14: lime evolution of the cloud water charge density at 

4 km TAGL ) following the lightning discharge. Solid lines represent 
positive charge while dashed lines represent negatively charged cloud 
droplets. 
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Fig. 15 : Same as in Fig. 14 except for the charge on graupel particles. 
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We know that the eneryy density produced by the electric field can 
be calculated from the following equation 


1 


- c E 2 
2 o 


(26) 


where e 0 = 8.86xl0" 12 F/m is the permittivity of air. If we integrate 
this expression over the whole domain of the simulated cloud using a 
numerical integration scheme, we can get the total electrical energy at 
each simulation time. We did this for Case 5, which represents a 
charge transfer of 4.8 Coul . A calculation over the cloud volume 
before the lightning discharge resulted in a total electrical energy of 
1.14xl0 10 J. The same calculation after the discharge gave a total 
energy of y.4xl0 9 J. Thus, the energy dissipated by the lightning 
channel is about 2xl0 9 J. If we consider the total simulated lightning 
channel length to be 428U meters (see Fig. 7) and assume that the 
energy dissipation along the channel is uniform, we calculate an energy 
dissipation along the channel of around 4.7x10 s J/m. 

Krider et al . (1968), based on observational data, inferred the 
energy dissipation along the lightning channel to be about 2.3xl0 5 J/m. 
Hill (1971) and Plooster (1971) developed models of the detailed 
lightning channel itself. In these models, a portion of the plasma 
channel was simulated without considering any of the environmental 
influences. From their modeling they calculated lightning energy 
dissipations of 1.5xlU 4 and 2.5xl0 3 J/m, respectively. There is an 
unresolved debate about the discrepancy between these two model calcu- 
lations and between the model calculations and the calculation made 
from laboratory simulation. Despite the simple numerical integration 
technique employed in our calculation and the assumption of constant 
energy dissipation along the channel, the model generated energy 
dissipation is in reasonable agreement with the work of Krider. 

Although the energy estimate does not agree so well with the works of 
Hill and Plooster, we find these initial estimates to be encouraging 
and add to our confidence in the parameterization scheme. 

4 . 4 Clou d -to-Ground Discharge Parameterization 

Having made these tests of the lightning parameterization scheme, 
we undertook a preliminary implementation of a cloud-to-ground 
discharge. Since our termination criterion of an electric field 
strength less than 160 kV/m would not allow a ground flash to occur 
in the 19 July simulation (and in fact no such ground discharge was 
observed by the lightning detection network installed for the project), 
we had to artificially force the intracloud discharge channel to go 
to ground by ignoring the lower termination criterion. The remaining 
initiation, propagation, and termination criteria were employed as 
for the intracloud discharge. The parameterization scheme has to be 
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modified somewhat, however, when the lightning channel makes 
contact with the ground, particularly with respect to the charge 
redi stri buti on. 


As our first order approximation, we assumed that Eqs. (20), 

(21), and (25) were all applicable because the physical concept of the 
cloud-to-yround discharge is the same as for the intracloud discharge. 
The key parameter in simulating a cloud-to-ground discharge is the 
electric potential of the lightning channel. Because the lightning 
channel is a conductor, we have where is the potential of 

the initiation point of the lightning channel. Because the earth is 
a good conductor, as soon as the channel attaches to the ground, both 
the' channel and earth should have the same value of potential. We 
choose the ground as the reference point (zero potential) for the 
potential field so the channel experiences a sudden change in potential 
when ground attachment occurs. The total charge transferred to the 
ground due to the cloud-to-ground discharge can be written as 


Vi - 



dn 


(27) 


where the integral is taken along the channel path. 

For this simulation, we chose k=1.5xl0“ 11 F/m (case 7) as the 
framework within which to test the cloud-to-ground discharge. Using 
Eq. (27), we obtained a transfer of -1.5 C of charge to the ground. 

This is within the observed range of charge transfer values for a 
cloud-to-ground discharge. Figure 16 shows the time evolution of the 
total charge density and cloud water charge density for a one minute 
period following the simulated discharge. Figure 17 shows the same 
history for the horizontal and vertical electric field components. 

As seen in Fig. 16, the lower part of the channel creates a large 
number of positive ions which attach to and reverse the polarity of a 
portion of the negative cloud water charge structure of the lower 
cloud. Initially, the total charge density field is dominated by the 
lightning produced ions, although over the course of the minute 
following the discharge, the original charge structure is seen to become 
more prevalent. Shortly after the simulated discharge, the original 
negative cloud charge nearly disappears and a new positive space charge 
is created on the cloud water. This positive charge accumulation 
strongly influences the vertical electric field component, as can be 
seen in Fig. 17. The positive vertical field region present beneath 
the cloud prior to the discharge disappears. A negative vertical field 
becomes dominant below the cloud after the lightning. The horizontal 
electric field component is also dramatically changed. 

Observations typically show a reversal of the electric field (as 
measured at the ground) from foul weather (positive) values to fair 
weather (negative) values following a lightning discharge. This abrupt 
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field reversal is followed by a gradual recovery to foul weather values 
if no other lightning occurs near in time to the discharge responsible 
for the initial reversal. As can be seen from the maximum negative 
values of the vertical field component in Fig. 17, a recovery process 
is underway in the negative field region below the cloud. While this 
gives qualitative agreement with an observed characteristic of 
cloud-to-ground lightning, the magnitudes predicted by the model may 
not be comparable to those observed. However, we must keep in mind 
that, in this simulation, we forced a ground discharge. If one had 
been predicted to occur naturally, the results may have been more 
meaningful . 

Another aspect of the parameteri zed discharge which appears in the 
simulation and has been observed to occur in actual thunderstorms is 
the creation or enhancement of the lower positive charge center 
following lightning. Holden et al . (1983) and Marshall and Winn (1982) 
reported on the creation of lower positive charge centers in 
thunderstorms in New Mexico following lightning discharges. They found 
that positive charge was deposited on precipitation falling from the 
base of the storms which subsequently fell to the ground. Their 
instrumentation was only capable of measuring the charge on 
precipitation-sized particles. In the model results shown in Fig. 16, 
we can clearly see the formation of this positive charge in the total 
charge density. Most of the charge is attached to the cloud particles 
(as noted above); however, there is also attachment of positive 
lightning ions to the rain falling in the vicinity of the discharge. 
This is further indication that the basic physical processes involved 
with the lightning discharge are being simulated in the model. The 
question of correct magnitudes remains to be investigated in more 
detai 1 . 

4.5 EMA Subcontract 


During the third year of this investigation, we engaged Electro 
Magnetic Applications, Inc., of Denver, CO, under a subcontract to aid 
in improving two of the criteria used in developing the lightning 
parameterization scheme. The four criteria involved in the scheme 
include initiation, propagation, termination, and charge redistribution 
which have been described in Section 4.2. The first of EMA's efforts 
involved an attempt to improve the propagation scheme for the 
discharge. The second effort involved investigating the effect of 
thunderstorm particles on the initiation criterion. 

4.5.1 The lightning propagation model 

Under the subcontract, EMA developed a Lightning Discharge 
Propagation Model (LDPM) for incorporation into the SEM. The purpose 
of the LDPM is to account for the tortuous path that is usually 
followed by a lightning channel. One of the weaknesses of the SEM 
discharge parameteri zation scheme has been the propagation criterion 
which artificially constrains the channel to follow the side of a grid 


- 49 - 



- 50 - 


ORIGINAL PAGE IS 
OE POOR QUALITY 


Fig. 16: Plots of the total charge density (A) and cloud water charge density (B) at various times 

following the cloud-to-ground discharge. Solid lines represent positive values and dashed lines, 
negative values. Maximum and minimum values are indicated on each plot. 
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Fig. 16 : (concluded.) 

















box, or its diagonal, depending on the ambient electric field direc- 
tion. The LDPM scheme takes into account the electric fields created 
by the leader tip and the previous leader segment. Also included are 
the effects of random free electrons in the vicinity of the leader tip 
as an influence in determining the direction of propagation of the 
next segment. 

The basic idea employed involves consideration of conditions in 
proximity to the leader tip. The leader tip contains a high 
concentration of charge. In the region around the tip, free electrons 
are randomly generated by the action of cosmic rays. These electrons 
are accelerated in the combined electric fields of the ambient 
thunderstorm and the leader. If the acceleration is strong enough, 
such electrons can collide with and ionize neutral air molecules 
creating secondary free electrons which are subsequently influenced by 
the field. As the process proceeds, an avalanche of electrons will 
occur and the air will become highly conductive in the avalanche 
region. The LDPM algorithm determines which free electrons will result 
in the conducting path that will become the next leader step. By using 
random seeds, the LDPM produces lightning channels which propagate in 
the same general direction as the original scheme used in the SEM, but 
show deviation from the path more reminiscent of actual lightning. The 
scheme has been included in the SEM but not as yet tested on an actual 
case. A report describing the LDPM in more detail is included in 
Appendix B. The report includes several figures which show the 
difference between the original propagation scheme and the improved 
scheme. 

4.5.2 Particle effects on lightning initiation 

The initiation criterion used for the SEM requires an electric 
field strength in excess of 400 kV/m in the cloud domain. In clear air 
at altitudes where we would expect to find lightning, the electric 
field required to initiate air breakdown would exceed 1000 kV/m. The 
model criterion of 400 kV/m is a rough estimate of what the actual 
breakdown field might be inside a thunderstorm based on observations of 
in-cloud electric fields that have rarely exceeded 400 kV/m. What the 
actual field strength is at the actual point of lightning initiation is 
an open question which may never be precisely answered. However, the 
evidence does suggest that the in-cloud breakdown field is considerably 
lower than the clear-air value. One of the most appealing explanations 
for this discrepancy is that the hydrometeors within a thunderstorm, 
particularly the ice particles, act to enhance the local thunderstorm 
field to the point where breakdown potential is achieved. This 
conjecture formed the basis for the EMA study of particle field 
enhancement . 

Their study involved an examination of how the shape, size, and 
number density of particles at different altitudes altered the pure air 
breakdown value. They accomplished this by looking at how the electron 
avalanche rate was affected by the various particle parameters. Since 
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analytical solutions for the electric field enhancement around 
irregularly shaped particles such as snowflakes and ice crystals do not 
exist, they employed an approximation using the field distribution of a 
spherical conductor (dipole field) and modified it with additional 
enhancement factors determined from previous research. Their results 
showed that only the largest uncharged ice particles were capable of 
making significant reductions in the breakdown field. When these 
particles were charged, the reduction of the breakdown field showed a 
linear dependence with particle charge and the enhanced field became 
asymmetric, favoring one-sided corona. While not studied per se, the 
interactions of two or more particles in further enhancing the local 
field were noted and suggested as being in need of further research. A 
report covering this work is also included in Appendix C. A scheme for 
incorporating these results into the initiation criterion of the SEM 
will be worked out at some future time. 
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5. SUMMARY AND CONCLUSIONS 


The results of this research have served two purposes. The first 
is gaining an insight into the low level strike environment for the 
F 106 Storm Hazards Project and the second is the development, for the 
first time, of a lightning parameterization scheme within the context 
of a multidimensional cloud electrification model. 

With regard to the low level strike environment, the simulation of 
the Wallops Island case has shown that, within the limitations of the 
model physics, a region of high electric field strength does occur in 
the 15-2U kft altitude range. In the early developmental stage of the 
storm, this high field region is associated with the presence of 
graupel , and possibly turbulence, making it a region where flight 
operations would have to be restricted. However, it is entirely 
possible that at a later time in the storm's life cycle, such a high 
field region would continue to exist without being accompanied by such 
threatening conditions. In addition, the Wallops Island simulation was 
characteri zed by little wind shear, leading to a nearly vertical 
development of the modeled storm. In a case where a stronger 
environmental shear were present, the high field region might not be 
concomitant with the presence of graupel or the graupel production 

might not be as heavy as in the case simulated. As a result of this 
investigation and complementary radar studies, an operational procedure 
based on these considerations was adopted for the 1986 field season. 

The development of the lightning parameterization scheme has been 
pioneering in nature. The simulation of a lightning channel within the 
context of a numerical cloud model has not been attempted before. We 
accomplished this by establishing criteria for the initiation, 
propagation, termination and charge redistribution of a lightning 
channel. Even though the parameterization is crude and the lack of 
understanding of how the lightning mechanism acts in a cloud is great, 
the simulations and studies we have done with the scheme active in the 
model have yielded results which are qualitatively consistent with the 
observations which are available. Because of the nature of this aspect 
of the research, more questions have been raised than have been 
answered; however, it appears that we have made a positive first step 
upon which to build as our understanding of lightning improves. Much 
more research into the storm-lightning environment is needed. With the 
lightning parameterization and the improvements developed through the 
EMA collaboration included in the SEM, it should be possible to make a 
significant contribution to the understanding of the thunderstorm 
envi ronment . 
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APPENDIX A 


THE PARAMETERIZATION OF THE LIGHTNING DISCHARGE IN A 
TWO-DIMENSIONAL, TIME-DEPENDENT ELECTRICAL CLOUD MODEL 

John H. Helsdon, Jr. 


ABSTRACT 


In order to extend the capabilities of thunderstorm electrification 
models to simulate all but the earliest stages of a storm, it is 
necessary that the electrical stresses created by an ever increasing 
electric field be relieved in some manner. In nature, it is the 
lightning discharge which fulfills this need. In keeping with the 
modeler's goal of attempting to simulate natural processes as faith- 
fully as possible, it is necessary for the model, then, to contain 
some parameteri zation of the lightning discharge. 

In this presentation, a first order attempt at simulating the 
intracloud lightning discharge within the context of a two-dimensional, 
time-dependent thunderstorm electrification model will be discussed. 

The criteria for the initiation, direction of propagation, termination, 
and redistribution of charge associated with the simulated channel will 
be presented. Finally, the effect of the parameterization in an actual 
simulation will be shown. 



A LIGHTNING PARAMETERIZATION SCHEME IN A TWO-DIMENSIONAL, 
TIME-DEPENDENT STORM ELECTRIFICATION MODEL 


John H. Helsdon, Jr., Richard D. Farley, 
and Gang Wu 


ABSTRACT 


An important consideration in the modeling of thunderstorm 
electrification is the existence of lightning in nature. A model 
without lightning is only capable of simulating the development of a 
storm in its early stages (up to the time of first lightning). Since 
lightning accounts for significant charge rearrangement within a 
storm, it becomes necessary to include such processes in any modeling 
effort which hopes to shed light on the electrical evolution of such 
storms. In an attempt to make allowance for the lightning discharge 
in our Storm Electrification Model (SEM), we have developed a 
simplified parameteri zation scheme to simulate the lightning channel 
and its attendant charge redistribution. 

This parameterization scheme has criteria for the initiation, 
propagation, termination, and charge redistribution associated with 
lightning. In order to simplify the scheme as much as possible, the 
first three criteria are all dependent on the strength and direction 
of the electric field vector. The propagation of the channel is 
calculated so as to maintain overall charge neutrality for an 
intracloud discharge. The charge deposited along the channel is then 
spread in a conservative manner laterally away from the channel over 
several grid points in order to maintain computational stability. 

This distribution of charges is then added to the small ion 
population and allowed to interact with the other charge carriers 
as the integration proceeds. 

A simulation of the 19 July 1981 CCOPE case has been run with the 
lightning discharge included. A parametric study has been undertaken 
to determine the effects of different amounts of charge transfer on 
the subsequent electrical evolution. The results of these simulations 
will be discussed. 



A NUMERICAL MODELING STUDY OF A MONTANA THUNDERSTORM: 
2. MODEL RESULTS VERSUS OBSERVATIONS 
INVOLVING ELECTRICAL ASPECTS 

John H. Helsdon, Jr., and Richard D. Farley 


ABSTRACT 


In this investigation, we compare the results of the Storm 
Electrification Model (SEM) simulation of the 19 July Cooperative 
Convective Precipitation Experiment (CCOPE) case study cloud against 
the actual observations with respect to the cloud's electrical charac- 
teristics, as deduced from the data of two aircraft. It is found that 
the SEM reproduces the basic charge and electric field structure of 
the cloud on a similar time scale to that observed. The comparability 
of the modeled and observed values of field strength and charge within 
the cloud is directly related to the proximity of the aircraft to the 
main active core of the cloud. The character of the electric field 
external to the cloud appears to be well modeled, at least at the time 
of observation. A feature which was not observed, but appears in the 
model, is a charge-screening layer at the top of the anvil. In addi- 
tion, the model predicts electric field strengths sufficient to 
initiate a lightning discharge at approximately the same time in 
cloud evolution as that when an intracloud discharge was recorded. 


A-4 



APPENDIX B 


EMA-86-R-55 


A STOCHASTIC MODEL FOR THE TORTUOUS 
PROPAGATION OF A LIGHTNING STEPPED LEADER 


by 

Dale A. Steffen 
Robert W. Haupt 
Terence Rudolph 


Prepared for: 

South Dakota School of Mines 
Institute of Atmospheric Sciences 
501 E. St. Joseph Street 
Rapid City, South Dakota 57701-3995 


Under Contract No. E02-86-001 


September 1986 


B-l 


A STOCHASTIC MODEL FOR THE TORTUOUS PROPAGATION 
OF A LIGHTNING STEPPED LEADER 


The results of Task 1 produced by EMA for the referenced contract 
"Completion of a Lightning Discharge Propagation Model (LDPM)" are discussed in 
this report. A FORTRAN code has been written for this model which can be integrated 
into the South Dakota School of Mines (SDSM) time-dependent Storm Electrification 
Model (SEM). The LDPM has been structured to enable its use in either a two or 
three-dimensional SEM with minor modifications. 

1. INTRODUCTION 

The LDPM is an approximation of the tortuous path of a lightning channel in 
a thunderstorm environment. Location and direction for the initial discharge point must 
be supplied as input data to the model. The discharge channel is then established as 
a series of connected stepped leader segments and propagates until certain 
termination conditions are met. All leader segments are assumed to be the same 
length and have equal charge densities. These parameters can be set to any 
values desired. While testing the LDPM, parameter values of 50 meters and 
1 millicoulomb/meter were used based on the average values published by Uman 1 . 

The primary parameter calculated in the model is the direction for a segment 
of the stepped leader. This computation is based on the effects of random free 
electrons in the vicinity of a leader segment tip which determine the direction of the 
next stepped leader segment. The details of the direction parameters for a stepped 
leader segment will be discussed in the next section. Once this direction is 
established, position vectors for a segment can be calculated. The model is then able 
to continue the computation of the positions and directions of following stepped leader 
segments. 


1 Uman, M.A., lightning. ’ McGraw-Hill, New York, 1969. 
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The remainder of this report discusses the calculations used in the LDPM, 
presents results obtained during the model tests, and provides a discussion on how 
the code is used. In addition, a listing of the code is provided. 

2. PHYSICAL CONSIDERATIONS 

Uman 1 discusses several theories which have been presented since the 
late 1 930's to explain observed stepped leader characteristics. Although basic 
differences exist between these theories, several fundamental issues are common for 
the theories and are listed below. 

1 . A stepped leader has an inner conducting core wherein most of the 
current flows. 

2. The core is surrounded by a corona region in which the electric field 
exceeds the breakdown field of air (approximately 3 x 10 6 volts/meter 
for standard temperature and pressure at sea level). 

3. A mechanism exists which precedes the stepped leader to create an 
ionized channel in "virgin air" which provides continuation of the 
stepped leader propagation. 

4. As current flow increases into the ionized channel, a condition occurs 
whereby the conductivity increases until the channel becomes a new 
segment of the stepped leader and the end of the channel becomes the 
new tip of the leader. 

These theories agree that stepped leader segment directions are randomly 
distributed. Based on lightning channel tortuosity studies by Hill 2 , the mean absolute 
direction change is nearly constant from flash-to-flash and on the order of 16 degrees. 
For purposes of the LDPM, it is assumed that the randomness of the stepped leader 
direction is related to the random creation of free electrons in air due to cosmic ray 
bombardment. 

2 Hill, R.D., 'Analysis of the Irregular Paths of Lightning Channels.* J. Geophvs. Res .. Vol. 73, 1968 
(pp. 1897-1906). 
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Issue 3 is of primary concern for this model. The stepped leader 
propagation theories often vary when describing this mechanism. The process 
occurring which results in the formation of the ionization channel ahead of the 
stepped leader tip has been more recently described by Rudolph 3 from studies 
involving the interaction of lightning with aircraft in flight. 

The leader tip is a region which contains a large charge density and high 
electric field levels. Free electrons (generated by cosmic rays) in the air are 
accelerated in this high field until they collide with a neutral atom or molecule. If the 
electron's kinetic energy is large enough at the time of the collision, the neutral particle 
can have an electron separated from it, producing a second free electron and a 
positive ion. The free electrons are then again accelerated by the field, possibly 
suffering more collisions and producing more free electrons and ions. If the rate of 
production of free electrons is larger than the rate of loss (by recombination, and 
attachment to form negative ions), an electron "avalanche" occurs, in which sufficient 
numbers of electrons and ions are produced and substantially alter the electrical 
conductivity of the air. 

Figures 1 and 2 illustrate how the effects of avalanching result in the 
generation of a new leader step from an existing leader tip. Figures 1 and 2 also 
illustrate the ionization path mechanisms for a negatively charged and for a positively 
charged leader tip, respectively. 

In Figure 1 , a free electron (at time ti) is accelerated by the electric field 
towards the leader tip. As avalanching occurs (time t 2 ), a net positive charge is formed 
in the region of the initial free electron position. This results from positive ions which 
have a much lower mobility or drift velocity than electrons. The net effect is 
neutralization of the electric field at the leader tip (time t 3 ). The conductivity in the 

region approaches the conductivity in the inner core of the stepped leader. Thus, the 
region in which positive ions were first produced becomes the new leader tip. 

In Figure 2, a free electron (at time t-j) is accelerated by the electric field 
away from the leader tip. As avalanching occurs (time t 2 ), a net positive charge is 

3 Rudolph, T., and R.A. Perala, "Linear and Non-Linear Interpretation of the Direct Strike Lightning 
Response of the NASA FI 068 Thunderstorm Research Aircraft,* Electro Magnetic Applications - Denver, 

EMA-83-R-21, March 21, 1983. 
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Figure 1. Time Scenario of Free Electron Avalanching for 
Positively Charged Leader Tip 
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Figure 2. Time Scenario of Free Electron Avalanching for 
Negatively Charged Leader Tip 


formed in the region of the initial position of the free electron. Moreover, as the 
electrons move away from this region, a conducting channel is generated between the 
positive ions and electrons until the electric field between the two opposite charge 
regions is neutralized. However, the electric field is enhanced between the leader tip 
and the positive charge region and induces additional forces on the negative charge 
region around the tip. In addition, the breakdown region is extended toward the 
positive charge region. This effect is illustrated in Figure 2 at time t 3 as a new 

conductive region extends the stepped leader. The leader tip occurs at the negative 
charge region created by the electrons which were first accelerated away from the 
original free electron position. 

The mechanism for extending the stepped leader is dependent on the 
location of free electrons which create a conductive channel to the leader tip. Thus, it 
is necessary to consider all free electrons and determine which one will initiate the first 
complete ionized path to the leader tip. Moreover, the algorithm used must determine 
which electron will initiate this complete ionized path in the minimum time 
(compared to other electrons). 

A rigorous mathematical treatment of the mechanisms just discussed would 
require involved calculations for the avalanche process and resultant conductivity in 
the stepped leader tip vicinity. These calculations would also involve significant 
computer time when used in the SEM. A simplified algorithm which has similar 
dependencies on the position of random electrons and the electric field around the tip 
of a stepped leader segment was developed. 

Since the drift velocities of the positive ions are significantly lower than those 
of electrons, the first simplification does not account for the mobility of resultant positive 
ions. The equation of motion for electrons generated by the avalanche is then given 
by: 


dV e 

dt 


Pe_ 

m e 


(E T + V e x B ) - 



where: q e 

m e 


electron charge, 
electron mass, 


(1) 
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Ej = total electric field, 

V e = electron velocity, 

B = magnetic field, and 
v = an approximate collision frequency. 

To further simplify this equation, it is assumed that the creation of electrons 
due to avalanching has progressed to a level where a collision dominated, plasma 
region exists. In turn, the electron velocity is constant and equal to the drift velocity, 
and: 


dV e 

dt 


= 0 . 


( 2 ) 



where: p e = electron mobility, and 

V<j = drift velocity. 

It is assumed that the vector r e between the leader tip and the electron 
position can be approximated by: 

r e = V d t, (4) 

or the time necessary to create a channel can be approximated: 

, . . iii . (5) 

I Vd I | Ej | 

It is assumed in the model that the ionized channel induced by the 
avalanching process (resulting from randomly occurring electrons) will initiate the next 
stepped leader segment in the same direction for either a positively or negatively 
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charged stepped leader segment. Thus, the model is independent of the leader tip 
charge. Based on these assumptions, the new stepped leader segment direction is 
determined by identifying the free electron in the vicinity of the previous stepped 
leader segment tip which will create an ionization path in a minimum time (compared 
to other free electrons). 

The model requires the following computations which will be discussed in 
detail in the next sections. 

a. Random generation of electrons in the vicinity of each stepped leader 
tip. 

b. Total electric field strength and direction at each electron location. 

c. The direction vector of the new stepped leader. 

3. LIGHTNING DISCHARGE PROPAGATION MODEL 

CALCULATIONS 

3.1 Description 

The LDPM code was written to predict the tortuous path of a lightning 
channel based on the physical considerations in Section 2. The routine is 
programmed to calculate the parameters for each stepped leader segment. 

3.1.1 Random Electron Generation 

A Monte Carlo technique was used to generate the random locations of free 
electrons within a sphere around each leader segment tip. The ambient electron 
density is assumed to be .2 electrons/m 3 . The size of the sphere is dependent upon 
the number of free electrons. Each model run in this report considered 500 random 
electrons around each tip. Thus, the sphere centered at the leader tip was required to 
have a radius of 8.419 meters. 
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A seed integer is also required for the Monte Carlo routine. This seed 
controls the random placement of electrons. Different seeds will produce different sets 
of electron locations. 

3.1.2 Total Electric Field Calculation 

The total electric field strength and direction was computed for each electron 
location generated by the Monte Carlo technique. The electric field at any electron 
location should be the resultant of the ambient field and the fields due to the previous 
leader segments. The model code was tested to incorporate the effects of the previous 
leader segments and the ambient field in the electric field calculation. It was observed 
that this calculation only required the ambient electric field and the electric field 
components due to the two leader segments closest to the electron. Electric fields 
caused by segments further away did not affect the modeled lightning path 
significantly. The total resultant electric field at a location (x.y.z) in the vicinity of the ith 
stepped leader segment is stated below: 

^total ( X >Y' Z ) = ^ambient ( x >y* z ) + segment ( x >y>^) + ^i-1 segment ( x »y>^) • (8) 

The general equation for the electric field at any point in space due to a line 
charge is given below and is illustrated by Figure 3. 

The electric field at any point in space due to the line charge element Xdl is: 


Xd\ (r - 0 

dE = ; X. = charge/unit length . 

4jte 0 I r - ? 1 3 

Thus, the electric field at point p due to the total line charge is: 

i = _l_ r dl . 

4n£ 0 Jo I r - (R + lu) | 3 


(7) 


(8) 
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Figure 3. Stepped Leader Segment Line Charge and 
Its Parameters in World Coordinates 


The general world coordinate solution for the electric field due to the line charge 
becomes: 


/ 

p -i 


p- _ 


X 

A(L - B) + u(A 2 - BL) 


A 2 u - AB 


4tc£ 0 (A 2 - B 2 ) j 

(A 2 - 2BL + L 2 ) 1 / 2 


IaI 


[ 




/ 


where: A = r - R , and 

B = A • u. 
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Figure 4 depicts the constant electric field contours due to a line charge. 
The arrows show the direction of the electric field. 


3.1.3 Direction Vector of the New Stepped Leader Segment 

Only one of the electrons generated in the vicinity of the leader tip will 
properly satisfy the model requirements and determine the direction of the new 
stepped leader segment. For simplification, each stepped leader segment is assumed 
to be positively charged (as discussed in Section 2) for the LDPM. The new stepped 
leader segment is assumed to be initiated by the ionized path formed in the least time. 
The time to generate the ionized path is proportional to the ratio formed from the 
distance of the leader tip to the electron (r e ) and the resultant electric field strength (Ej) 

at the electron location. This relationship was determined in Section 2 and is stated in 
Equation 5. Thus, the new stepped leader segment direction is chosen to be the 
electric field direction at the initial location of the free electron which initiates an 
ionized path in the minimum time (compared to other electrons). 


If the angle 0 e , formed by the Ej and r e vectors at the initial electron location 
is within a certain angle range, that electron will create an ionized channel to the 
leader tip along a "relatively" straight path. The time it takes to establish the ionization 
channel can be approximated using Equation (5). As 0 e increases beyond this certain 

angle range, the actual path of the ionization channel along electric field lines will 
deviate from a straight line. Therefore, the actual time it takes to establish the ionized 
channel to the leader tip may be significantly greater than the time approximated by 
Equation (5). In fact, if 0 e is large enough, the ionized path may not attach to the 
leader segment tip, but to some point behind the tip instead. The cosine of 0 e formed 
by Ej and r e at any electron location is computed as follows: 


cos 0 e 


* E T 

r e l IE t I 


( 10 ) 


A stepped leader segment is defined by the computer code as a finite length 
with its tip being a single point. In nature, the leader tip is not a single point but is 
actually defined as a volume. In terms of the LDPM, this volume exists around the 
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Figure 4. Constant Electric Field Amplitude Contours for Line Charge 
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leader tip point and has an electric field strength greater than 3 x 10 6 V/m everywhere 
inside. 


Figure 5 depicts the above concepts. The angle range parameter is 
determined by the radius of the corona region. This radius is on the order of 2 meters 
and is assumed to be the same magnitude for all leader segments in the LDPM. The 
angle range parameter was adjusted to an absolute value of 12.5 degrees. This value 
appeared to allow the lightning channel path to be reasonably tortuous. The angle 
range parameter can be changed in the LDPM as desired. However, 12.5 degrees is 
recommended, since smaller angle parameter values limit the tortuousity, while larger 
angle parameter values result in erratic tortuousity which exceeds the findings of Hill 2 . 

Each of the 500 electrons generated by the Monte Carlo routine is tested in 
the code to determine which one is used to establish the direction of the new stepped 
leader segment. The one electron chosen must meet the following conditions: 

1 . The initial free electron which initiates an avalanche must originate in a 
location where the electric field strength is less than 3 x 10 6 V/m. (If the 
field was larger, the electron was assumed to be inside the corona 
region of the leader segment and was Qfii considered). 

2. The angle 0 e formed by Ej and r e for an electron location must be less 

than the angle range parameter (absolute value of 12.5 degrees) 

(If 0 e was greater, the electron was rigl considered). 

3. Of the remaining electrons, the electron which creates an ionization 
channel to the leader tip in the least time, t (as computed in 
Equation (5)) was identified and used to calculate the new direction. 

The new stepped leader segment direction is then taken to be the direction of the total 
electric field at this "identified" electron location. 

It is important to note that, if the Monte Carlo sphere is too small, most or all 
of the random electrons will lie in a region where the electric field strength is above 
3 x 1 0 6 V/m. In turn, the model may not function properly. It was observed that a 
sphere containing 500 electrons and having a radius of 8.419 meters worked best 
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Electron "C" will not be drawn to the leader tip at all and equation [5] does not apply. 

However, e a >12.5* and the electron will not be considered. 

Figure 5. Paths of Free Electrons in Vicinity of Stepped Leader Tip 
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during the LDPM development. This sphere size contained enough electrons in a 
region where the electric field was below breakdown and allowed the LDPM to 
function properly and efficiently. 

3.2 Model Termination Conditions 

The LDPM will terminate if any of the following conditions are met: 

1 . The stepped leader enters a defined region (grid cell) with an ambient 
electric field less than 1.5 x 10 5 V/m. This termination criterion was 
based on work done by Griffiths 4 involving the propagation of positive 
streamers in the laboratory, 

2. The stepped leader propagates outside the designated problem space, 

3. (Optional) The stepped leader enters a defined charged region 
(Specified grid cell). 

3.3 MODEL INPUT AND OUTPUT REQUIREMENTS 

3.3.1 Input 

The problem space is depicted as a rectangular grid either in two or three 
dimensions. The code is written to accommodate two or three-dimensional data. This 
version of the LDPM is adjusted for a two-dimensional SEM. Each grid cell has side 
lengths equal to 200 meters. This routine requires the number of grid cells along the 
x-axis (horizontal), along the z-axis (vertical), and the initial discharge coordinate (x,z). 
Ambient electric field components are required for each grid cell. The routine also 
requires a seed integer for the Monte Carlo scheme. These variables are: 

1 . NCELX Number of cells along the x-axis 

NCELZ Number of cells along the z-axis, 


4 Griffiths, R.F. and C.T. Phelps, "The Effects of Air Pressure and Water Vapor Content on the Propagation 
of Positive Corona Streamers and Their Implications to Lightning Initiation,* Quart. J. Rov. Meteor. Soc .. Vol 

102,1976 (p. 419). 
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2. 

ICX, ICZ 

Initial discharge point (x,z), 

3. 

FEX(K), FEZ(K) 

One-dimensional array of x,z ambient field 



components, and 

4. 

ISEED 

Monte Carlo seed integer; must be less 



than 32,000. 


For the NCAR Electricity Model: 



1 . 

NCELX = 97; 

NCELZ = 97 


2. 

ICX = 86; 

ICZ = 31 


3. 

FEX(K),K = 1,9409; 

FEZ(K),K = 1,9409 

3.3.2 

Output 



The routine generates the (x,z) position in world coordinates for each 
stepped leader segment tip in the propagation path. The x and z arguments are 
written in the one-dimensional RX and RZ arrays. Coordinate values x and z are given 
in meters from the origin (x = 0, z = 0). These arguments can be passed back into the 
main program and are as follows: 

1 . RX(NL), RZ(NL) One-dimensional arrays of the x and z values in 

meters 

NL corresponds to the leader number (not the 
grid cell number), and 

2. NFL Integer: total number of NL leader segments. 

3.4 MODEL EXAMPLES 

Several lightning discharge propagation paths were generated and are 
shown in Figures 6 through 8. Figure 6 shows the lightning paths modeled for the 
NCAR electricity data at 58.5 minutes. The solid and dashed lines are the models 
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Figure 6. NCAR Electricity Model Data (58.5 Minutes) 
(Comparison of Monte Carlo Seeds) 
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Figure 7. EMA Electricity Model Data - Initial Discharge 
in the Positively Charged Grid Cell 
(Comparison of Monte Carlo Seeds) 
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Figure 8. EMA Electricity Model Data - Initial Discharge 
in the Positively Charged Grid Ceil 
(Comparison of Different Ambient Electric Fields) 
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produced from the LDPM and the algorithm described in Helsdon 5 , respectively. The 
four plots in Figure 6 demonstrate the effects of varying ISEED (Monte Carlo Seed for 
Generating Random Electron Locations). 

Figure 7 shows the effects of ISEED for ambient electric fields generated 
from an in-house routine. This routine creates electric fields from a set of discrete 
charge centers. Figure 8 shows the effects of varying the charge magnitudes at 
specified cells. These ambient fields were computed from the in-house routine. 

4. MODEL VERSIONS AND DATA 

The listings for three versions of LDPM are included in Appendix A of this 
report. In addition, the code used to generate ambient electric fields is included. The 
listings are as follows: 

1. LDPM - Subroutine 

Subroutine with common statements; requires the number of x cells 
(NCELX), z cells (NCEL2), the initial discharge coordinate, ambient 
electric field components at each cell, and ISEED. Passes leader tip 
locations RX(I), RZ(I) and NFL, the number of total leaders back into 
main program. 

2. LDPM.FOR 

Routine reads in ambient electric fields; requires NCELX, NCELZ, 
initial discharge coordinate, and ISEED. Calculates leader tip 
locations. 

3. LDPM - Interactive (not required by contract; used as test program) 
Interactive routine; generates ambient electric field components over 
grid from specified cells and their charge magnitudes. Requires 
NCELX, NCELZ, the initial discharge coordinate, and ISEED. Routine 
also requires coordinates of specified charge cells. Calculates leader 
tip locations. 

5 Helsdon, J.H., R.D. Farley and G. Wu, "A Lightning Parameterization Scheme in a Two- 
Dimensional, Time-Dependent Storm Electrification Model," presented at the 1986 Conference 
on Cloud Physics, held Sept. 22-26, 1986 in Snowmass, Colorado. 
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4. EFS.DT 

NCAR ambient electric field data (58.5 minutes, unformatted) 
FEX(K), K = 1 , 9409; FEZ(K), K = 1 , 9409. 
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SUBROUTINE LDPM 

C DEVELOPED BY E.M.A. - AUGUST 1986 

C SUB-PROGRAM DETERMINES DIRECTION OF NEXT LEADER SEGMENT DUE TO 
C ELECTRON WHICH WILL FORM AN IONIZED PATH TO 
C THE STEPPED LEADER SEGMENT TIP IN THE LEAST TIME 

E FIELD TOTAL = E FIELD LEADER + E FIELD PREVIOUS LEADER 
+ AMBIENT E FIELD 

NCEL- • CELLS IN X AND IN Z 
ICX#ICZ= RXfRZ COORD OF INITIAL CHARGE REGION 
NFL= TOTAL NUMBER OF STEPPED LEADER SEGMENTS 

CQMMQN/MAIN/fiX(0:500)#RZ'0:500>#FEX(9409>»FEZ(?40?> 
COMMON/VAR/NCELX # HCELZ # ICX» ICZrNFL 
CHARACTERS OFILEl»OFILE2 

ISEED=RANDOM FUNCTION INTEGER 
ISEED=30997 

OFILE1 3 'LDPM.DAT' » BINARY PLOTFILE OF LIGHTNING PATH 
0FILE2= 'LDPM .PRINT' ! PRINT VALUES OF LIGHTNING MODEL 

INPUT VARIABLES 

CL=SEGMENT LENGTH ELAMBDA=LEADER CHARGE DENSITY 
BDF=BREAK DOWN FIELD EK 3 1/4PIEpso NL*LEADER ♦ 

DZ=Z CELL HEIGHT DX=X CELL WIDTH 

NE=# FREE ELECTRONS ED-ELECTRON DENSITY U2/®**3) 

Q=hAX ALLOWABLE ANGLE BTU EFIELD l DIRECTION VECTORS 

DATA CL # ELAMBDA r BDF/50 . t * 001 # 1 *5E5/ 

DATA PI #EKrNL/3. 1415927# 9 *E9#0/ 

DATA DZ>DXfNE»EDfQ/200. ? 200. #500* .2# 12*5/ 

DATA EXM1#EYM1iEZM1/3*0./ 

TUQPI 3 2.*PI 
GP=Q*P 1/180 * 

CONST =EKtELAMBDA 
VOL=FLOAT(NE)/ED 
RADIUS 3 < . 75WJL/PI >**< 1 ./3. ) 

KTOT=NCELX*NCELZ 
XTOT=NCELX*DX 
C 

C INITIALIZE STEPPED LEADER PROPAGATION 
C INPUT DATA: Rif AMBIENT E FIELD# U1 
C COORDINATES OF INITIAL POINT ON GRID PROBLEM SPACE 
RX(0)=DXt( FLOAT ( ICX)-*5) 

RY=0. ! FOR 2D CASE# Rr<0)*DY<FL0AT<ICY)-*5>r FOR 3D CASE 

RZ (G)=DZ1 (FLOAT ( ICZ) *.5) 

C 

IC=ICX+NCELX*( ICZ-1 ) 

C FOR 2D CASE 

E TOT-SORT (FEX(IC)*FEX(IC)+FEZ(IC)tFEZ(IC) ) 

C FOR 3D CASE 

C ETOT=SQRT(FEX(IC)tFEX(IC)+FEY(IC)*FEY(IC)+FEZ(IC>*FEZ(IC)> 

C DEFINE UNIT VECTOR 
UX=FEX( IO/ETOT 

UY*0. f FOR 2D CASE* UY 3 FEY(IC)/ETOT# FOR 3D CASE 
UZ-FEZ( IO/ETOT 
C 

0PEN(1 # FILE=0F I LE 1 * FORM* 'UNFORMATTED' # STATUS® 'NEW' ) 

OPEN< 2 #FILE 3 QFILE2* STATUS 3 'NEW' ) 

C 
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C START LOOP 
2000 NL=NL+1 
C 

C DETERMINE CELL POSITION OF NEW LEADER SEGMENT TIP ON OVERLAY 

C USE APPROPRIATE AMBIENT E FIELD MAGNITUDE AND DIRECTION 

C K IS THE CELL NUMBER OF THE LEADER SEGMENT TIP LOCATION 
KX=>ANINT((RX(NL-l)/DX) + .5) 

KZ=ANINT < (RZCNL-1 )/DZ)+.5) 

K=KXtNCELX*!KZ-l) ! NEED KT FOR 3D CASE 
C FOR 2D CASE 

EF AMB=S0RT ! FEX ( K ) *FEX ( K ) +FEZ ( K ) »FEZ (K) ) 

FOR 3D CASE 

EF AMB=SQRT < FEX < K > *FEX < K ) +FEY C K) *FEY ( IO+FEZ < K) «FEZ < K) ) 
TERMINATE IF 

LEADER SEGMENT HAS PROPAGATED OUTSIDE THE GRID 
LEADER SEGMENT ENTERING CELL WITH AN AMBIENT FIELD LT l.SM V/M 
IF!K.LT.1.0R.K.GT.KTOT.OR. 

+ RX(NL-l) .LT.O.OR.RX(NL-l) .GT .XTOT)GOTO 3000 

URITE<2.») 'LEADER SEGMENT NUMBER' »NL 
URITE(2»*> 'CELL NUMBER '.K.' <X»Z) COORD'.KX.KZ 

WRITE(2.*)'T0T E FIELD '.ETOT.' AMB E FIELD ' .EFAMB 
URITE(2.t)'AMB EX '.FEX(K>»' AMB EZ '.FEZOO 

URITE!2.*)'X COORD '»RX(NL-1)»' Z COORD 'iRZ(NL-l) 

URITEI2**) ' 

URITE(1)RX<NL-1).RZ(NL-1) 

IFIEFAMB.LT. BDFIGOTO 3000 

SHORTEST TIME OF ELECTRON TO FORM IONIZED CHANNEL 
THIS TIME DETERMINES DIRECTION OF NEU LEADER SEGMENT 
TMIN=1.E20 

USE MONTE CARLO METHOD TO GENERATE RANDOM FREE ELECTRONS 
TO CALCULATE ELECTRON DISTRIBUTION NEAR LEADER SEGMENT TIP 
RXP=RX(NL-1 J+UXICL 

RYP=RY+UY*CL ! USE RY(NL-l) FOR 3D CASE 
RZP=RZ(NL-1)FUZ*CL 

DO 100 1=1. NE 

IN WORLD COORDINATES 

CALL RANDQfU ISEED.RAND) 

RE=RAND*RADIUS 

C RE=SGRT ( RANDtRADIUStRADIUS ) ! FOR ACTUAL RANDOM DISTRIBUTION 

CALL RANDOM! ISEED.RAND) 

TH=RAND*PI 

C TH=AC0S((RAND*2)-1) ! FOR ACTUAL RANDOM DISTRIBUTION 

CALL RANDOM! ISEED.RAND) 

PH=RAND*TUOPI 

RMX=RE*SIN(TH)*COS(PH) 

RMY=RE*SIN!TH)tSIN!PH) 

RMZ=RE«COS(TH) 

RM=SQRT < RMX*RMX+RMY*RMY+RMZ*RMZ ) 

X=RMX+RXP 

Y=RMYFRYP 

Z=RMZtRZP 

C 

C CGMPUTE E FIELD COMPONENTS FOR Ith ELECTRON 
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AX=X-RX<NL-1) 

AY=Y-RY ! USE RY(NL-l) FOR 3D CASE 
AZ=Z-RZ(NL-l) 

A=SQRT ( AX* AX4AY*AY4AZt AZ ) 

IFC A.LT. 1 .E-£)GQTO 100 

R=AX*UX4AY*UY4AZ*UZ 

AS 3 A*A 

BS=B*B 

IFCAS-BS.LT * 1 .E-6)G0T0 100 
P=CONST/(AS-BS) 

D=CL*CL-2.*B*CL4AS 
IF(B.LE.O. )GOTO 100 
D-SQRT(D) 

C 

C E FIELD COMPONENTS DUE TO LEADER SEGMENT AND AMBIENT FIELD 
EX=P*( < < AX*(CL-B)4UXS< AS-B*CL) )/D)- 
4 ( <UX*AS-AX*B)/A) )4FEXCK> 

EY=Pt( ( (AY*(CL-R)4UY*<AS“B*CL) )/D)- 
4 ( <UY*AS-AY*B)/A) ) »4FEY(K) FOR 3D CASE 

EZ=P*( < (AZ*<CL-B)4UZ*<AS-B*CL> >/D>- 
4 < <UZ*AS-AZ*B)/A) >4FEZ(K) 

C 

IFCNL.LE.DGOTO 1000 
AXMl=X-RX(NL-2> 

AYM1=Y-RYM1 ! FOR 2D CASE? RYMl=RY(NL-2> » FOR 3D CASE 
AZM1=Z-RZ(NL“2) 

AN 1 =SQRT C AXM1 *AXM1 4AYM 1 *AYH1 +AZH1 *AZM1 ) 

IF(AM1.LT.1.E-A>G0T0 100 

BM1=AXM1 *UXM14AYM1*UYM14AZM1*UZM1 

ASM1=AM1*AM1 

BSMl=BHltBMl 

IF(ASM1-BSM1 .LT. 1 .E-6)GOTO 100 
PM1 =CQNST / ( ASM1 -BSM1 ) 

DM1 =CL*CL-2 . *BMUCL4ASM1 
IF (DM1 .LE.O. )GOTO 100 
DM1=SQRT(DM1 ) 

C 

C E FIELD COMPONENTS DUE TO PREVIOUS LEADER SEGMENT 

EXM1=PM1*<((AXM1*(CL-BM1>4UXM1*<ASM1-BM»CL>>/DM1>- 
4* ( <UXM1*ASM1-AXM1*BM1)/AM1 ) ) 

EYM1=PM1*( ( ( AYM1 t (CL-BM1 )4UYNlt(ASMl-BMl*CL) )/DMl )- 
4 <<UYM1*ASM1-AYN1*BM1)/AM1>> 

EZM1=PH1*( ( ( AZMI t ( CL-BM1 )4UZM1*(ASN1-BM1*CL) >/DMD- 
4 ( (UZMlIASMl-AZMlfBMl )/AMl ) ) 

C 

C TOTAL E FIELD COMPONENTS FOR Ith ELECTRON 
1000 EX*EX4EXM1 
EY*EY4EYM1 
EZ=EZ4EZM1 

ET =SQRT ( EX*EX4EY*EY4EZ*EZ ) 4 • 01 
DL * < RMX*EX4RM Y*E Y4RMZ*EZ ) / ( RN*ET ) 

C 

C CHOOSE APPROPRIATE ELECTRON TO FORM IONIZED CHANNEL 
C TO THE LEADER SEGMENT TIP IN SHORTEST TIME 
IFCET.LT . (3.0E6) .AND. (RM/ET) .LT.TMIN. 

4 AND . ACOS ( ABS ( DL ) ) . LE . QP ) THEN 
TMIN=RM/ET 
ETOT^ET 
DIRX = EX 
D1RY=EY 
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C CALCULATE DIRECTION FOR NEXT LEADER SEGMENT 
DRX=CL*UX 
DRY=CL*UY 
DRZ=CL*UZ 
RYM1=RY 

RY=KY+DRY » FOR 2D CASEi RY < NL) *RY(NL-1 ) +DRY f FOR 3D CASE 
RX(NL )=RX(NL“1 ) +DRX 
RZ(NL)=RZ(NL-1 )+DRZ 
C 

C NEW UNIT VECTOR 

C UNIT VECTOR IS PROJECTED ONTO X*Z PLANE FOR 2D SPACE 
C FOR 3D CASE 

C DIR=SQRT(DIRX*DIRX+DIRY*DIRY+DIRZ*DIRZ) 

C FOR 2D CASE 

DIR=S0RT ( DIRX*DIRXFDIRZ*DIRZ) 

IF ( DIR . GE . 1 . E-6 ) THEN 

UXhl=UX 

UYM1=UY 

UZhl^UZ 

UX=DIRX/DIR 

UY=0» ! FOR 2D CASE! UY=DIRY/DIR> FOR 3D CASE 

UZ=DIRZ/DIR 
END IF 
C 

GO TO 2000 
C 

3000 NFL=NL 

IF(K.LT» 1 .OR.K.GT *KTOT .OR* 

+ RX(NL-1).LT.0.0R.RX<NL-1).GT*XT0T)THEN 
URITE<2>*> 'LEADER OUT OF BOUNDS' 

END IF 

IF(EFAMB .LT »BDF)THEN 

URITE(2f <) 'LEADER SEGMENT ENTERED REGION LT BREAK DOWN' 

END IF 
C 

CLOSE(l) 

CLOSE (2) 

RETURN 

END 

C 

C 

SUBROUTINE RANDOM < ISEED* RAND) 

RANDOM NUMBER GENERATOR 

DOUBLE PRECISION E31rE32fSIfSJ 
£ 31 * 2 . 1*31 
E32=2.**32 
SI=DFLQAT ( ISEED) 

SJ=MOD( 69069. *SI+1 • »E32) 

IF ( ABS(SJ) .LT.E31 )THEN 
SI=SJ 

ELSEIF(SJ.GE.E31 )THEN 
SI=SJ-E32 
ELSE 

SI=SJ+E32 
END IF 


DIRZ-EZ 
END IF 

100 CONTINUE 
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IF < SI .GE.O. )THEN 

RAND=SI/E32 

ELSE 

RAND=1 .+SI/E32 
END IF 

ISEED*INT(SI) 

RETURN 

END 


LDPM.FOR 
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C DEVELOPED BY E.M.A. - AUGUST 1986 

C PROGRAM DETERMINES DIRECTION OF NEXT LEADER SEGMENT DUE TO 

C ELECTRON UHICH UILL FORM AN IONIZED PATH TO THE SEGMENT TIP IN 

C THE LEAST TIME 

C E FIELD TOTAL = E FIELD SEGMENT + E FIELD PREVIOUS SEGMENT 
+ AMBIENT E FIELD 

C NCEL= # CELLS IN < AND IN Z 
C ICX.ICZ 3 X.Y COORD OF INITIAL CHARGE REGION 

PARAMETER t NCcLX=97»NCELZ=97 * ICX=86 » ICZ=31 ) 

DIMENSION FEX ( 9409 ) f FEZ < 9409 ) 

CHARACTER*40 INFILE.OFILE1.0FILE3 
C 

C ISEED=RANPOH FUNCTION INTEGER 
ISEED=208?7 

C SPECIFIED CHARGED REGION CELLS (OTHER THAN INITIAL CHARGED CELL) * 

DATA NCI .NC2/195. 18/ * 

C 

INF1LE='EFS.DT' • INPUT AMBIENT E FIELDS 

OF ILE1= 'LDPN3.DAT' ! BINARY PLOTFILE OF LIGHTNING PATH 

0FILE3='LDPM.PRINT' ! PRINT VALUES OF LIGHTNING MODEL 
C 

C INPUT VARIABLES 

C CL=SEGHENT LENGTH ELAMBDA= SEGMENT CHARGE DENSITY 

C 8PF=6REAK DOUN FIELD EK=l/4PIE»»so NL “LEADER SEGMENT * 

C DZ-Z CELL HEIGHT DX=X CELL WIDTH 

C NE=* FREE ELECTRONS ED=ELECTRON DENSITY <.2/»M3> 

C Q=MAX ALLOUABLE ANGLE BTU EFIELD X DIRECTION VECTORS 
C 

DATA CL.ELAMBDA.BDF/50.. .001.1.5E5/ 

DATA PI. EK.NL/3. 1415927.9. E9.0/ 

DATA DZ . DX . NE . ED . 0/200 . . 200 . . 500. .2 . 12 .5/ 

DATA EXM1.EYM1.EZM1/3W./ 

TU0PI»2.*PI 

QP=Q*PI/180. 

CONST=EK*ELAMBDA 
VOL=FLOAT(NE)/ED 
RADIUS 3 ( . 75HV0L/PI > ** < 1 ./3 . > 

KTOT=NCELX*NCELZ 

XTOT=NCELX*DX 

C 

C READ IN AMBIENT E FIELDS 

0PEN12. FILE-INFILE.FORM^'UNFORMATTED' .STATUS*'OLD' ) 

READ(2)(FEX(K).K-1.KT0T).(FEZ(K)»K“1.KT0T) 

CLOSE! 2) 

C 

C INITIALIZE LEADER SEGMENT PROPAGATION 
C INPUT DATA: Rl. AMBIENT E FIELD. U1 
C COORDINATES OF INITIAL POINT ON GRID PROBLEM SPACE 
RX=DXt(FL0AT(ICX)-.5) 

RY*0. 

RZ=DZ*(FL0AT(ICZ)-.5) 

C 

IC= ICX+NCELXt ( ICZ-1 ) 

ETOT=SQRT (FEX < I C > M2+FEZ ( 10 ** 2 ) 

C DEFINE UNIT VECTOR 
UX“FEX(IC)/ETOT 
UY=0. 

UZ=FEZ( IO/ETOT 
C 
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OPEN ( 1 , FILE*QFILE1 , FORM* 'UNFORMATTED ' t STATUS* ' NEW ' ) 

OPEN ( 3 1 FILE-0FILE3 » STATUS* ' NEW ' > 

C 

C START LOOP 
2000 NL*NL+1 
C 

C DETERMINE CELL POSITION OF NEW LEADER SEGMENT TIP ON OVERLAY 

C USE APPROPRIATE AMBIENT E FIELD MAGNITUDE AND DIRECTION 

C K IS THE CELL NUMBER OF THE LEADER SEGMENT TIP LOCATION 
KX*ANINT< CRX/DXH.5) 

KZ=ANINT ( (RZ/DZH.5) 

K=KX+NCELX*<KZ-1> 

EF AMB*SQRT ( FEX ( K ) M2+FEZ CK)*t2) 

C 

C TERMINATE IF 

C LEADER SEGMENT HAS PROPAGATED OUTSIDE THE GRID 
C LEADER SEGMENT ENTERING CELL WITH AN AMBIENT FIELD LT 1»5M V/M 
C LEADER SEGMENT ENTERS INTO CHARGED REGION 
IFCK .LT» 1 .OR .K.GT .KTOT .OR. 
f RX *LT .O.QR.RX.GT*XTQT)60TO 3000 
C 

WR I TE ( 3 » * ) NL » K > KX » KZ i RX » RZ f ETOT t EFAMB # FEX C K ) ? FEZ ( K ) 

PRINT*, NL,K,KX,KZ,RX,RZ,ETOTrEFAMB,FEX(K),FEZ(K) 

URITEC 1 )RX,RZ 
C 

IFCEFAHB.LT *BDF) GOTO 3000 

c 

C SHORTEST TIME PATH OF ELECTRON TO LEADER SEGMENT TIP 
C THIS TIME PATH DETERMINES DIRECTION OF NEW LEADER SEGMENT 
TMIN*l.E20 
C 

C USE MONTE CARLO METHOD FOR FREE ELECTRONS 
C TO CALCULATE ELECTRON DISTRIBUTION NEAR LEADER SEGMENT TIP 
RXP*RX+UX*CL 
RYP*RYH1Y*CL 
RZP*RZ+UZ*CL 
C 

DO 100 1*1, NE 
C 

C IN WORLD COORDINATES 

CALL RANDOMC ISEEDtRAND) 

RE=RAND*RADIUS 

C RE=SGRT{RAND*CRADIUS**2> ) 

CALL RANDOMC ISEEDtRAND) 

TH=RAND*PI 

C TH- ACQS ( RANDB2- 1 ) 

CALL RANDOHC ISEEDtRAND) 

PH=RAND*TWOPI 
RE*RAN C ISEED ) tRADIUS 
TH=RANC ISEED)*PI 
PH*RANCISEED)*2.tPI 
RMX=RE*SINCTH)*COSCPH) 

RMY=RE*SIN( THXSINCPH) 

RMZ=RE*COSCTH) 

X*RMX hRXP 
Y=RMY+RYP 
Z=RMZ+RZP 

RM=SQRT ( RMX«2+RMYM2+RMZ W2 ) 

C 
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C COhPUTE E FIELD COMPONENTS FOR Ith ELECTRON 
AX*X-RX 
AY=Y-RY 

az*z-rz 

assort < ax«2+ay**2+ az**2 > 

IF( A.LT. 1 .E-6)G0T0 100 

B*AX*UX+AY*UY+AZ*UZ 

AS=A**2 

BS a B**2 

IF<AS-BS.LT.1.E-6>60T0 100 
P-CONST/(AS-BS) 

D=CL**2-2.*B*CL+AS 
IF(D.LE.O.)GOTO 100 
D=SC1RT(D) 

C 

C E FIELD COMPONENTS DUE TO LEADER SEGMENT AND AMBIENT FIELD 
EX=Pt ( ( ( AX* ( CL-B ) +UX* ( AS-B*CL ) ) /D ) - 
+ ( (UX*AS-AX*B)/A) ) +FEX <K ) 

EY=P* ( ( ( A Y* ( CL-B ) +UY* ( AS-B*CL ) ) /D ) - 
+ ( (UY*AS-AY*B)/A) ) !+FEY(K) 

EZ=P*( ( ( AZ * ( CL- B HUZ* < AS-BJCL ) >/D»- 
+ ((UZ*AS-AZ*B)/A))+FEZ(K) 

C 

IFfNL.LE.l )GOTO 1000 
AXM1=X-RXH1 
AYM1=Y-RYM1 
AZM1-Z-RZN1 

ANl=StlRT<AXMl**2+AYMlM2+AZMl«2) 

IFIAMl.LT. 1.E-6IG0T0 100 

BMl=AXNl*UXMl+AYHl*UYHltAZMl*UZMl 

ASM1*AM1M2 

BSM1=BM1**2 

IF ( ASH1 -BSM1 . LT . 1 .E-6 ) GOTO 100 
PMl=CONST/( ASM1-BSM1 ) 

DM1=CL**2-2.*BM1*CL+ASN1 
IF(DMl .LE.O. )GOTO 100 
DM1=SQRT(DM1 ) 

C 

C E FIELD COMPONENTS DUE TO PREVIOUS LEADER SEGMENT 

EXM1=PM1*(((AXM1*(CL-BM1)*UXM1*(ASM1-BM1*CL))/DH1)- 
+ <(UXH1*ASN1-AXH1*BM1)/AM1>> 

EYM1=PM1*(((AYM1*(CL-BMI>+UYM1*(ASM1-BM1*CL))/DM1>- 
+ ( (UYM1*ASM1-AYM1*BM1)/AM1) ) 

EZH1=PM1*( ( ( AZhl*(CL-BMl )+UZHl*(ASMl-BNl*CL) )/DMl)- 
♦ ( (UZM1*ASM1-AZM1*BM1 )/AMl ) ) 

C 

C TOTAL E FIELD COMPONENTS FOR Ith ELECTRON 
1000 EX=EX+EXM1 
EY=EY+EYM1 
EZ=EZFEZM1 

ET=SQRT(EX**2+EY«2+EZ**2)+.01 

DL=(RMX*EX+RMY*EY+RMZ*EZ>/<RH*ET) 

C 

C CHOOSE APPROPRIATE ELECTRON TO BE DRAWN TOWARD 

C LEADER SEGMENT TIP IN SHORTEST TIME 

IFIET .LT. (3.0E6) .AND. (RM/ET) .LT. THIN. 

+ AND.ACOS(ABS'DL)).LE.QP)THEN 
THIN=RM/ET 
ETOT=ET 
DIRX=EX 
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DIRY=EY 
DIRZ=EZ 
END IF 

100 CONTINUE 
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c 

C CALCULATE DIRECTION FOR NEXT LEADER SEGMENT 
DRX=CL*UX 
DRY=CL*UY 
DRZ a CL*UZ 
RXM1=RX 
RYM1=RY 
RZMZ^RZ 
RX^RXFDRX 


RY=RY+DRY 

RZ=RZ+DRZ 

C 

C NEW UNIT VECTOR 

C DIR=SQRT(DIRX**2+DIRYM2+DIRZW2) 

DIR=SQRT<DIRXt*2+DIRZ**2> ! PROJ'N ONTO 2D GRID 
IF( DIR*LT* 1 *E-6)G0T0 2000 
UXM1=UX 
UYM1=UY 
UZM1=UZ 
UX=DIRX/DIR 
C UY=DIRY/DIR 

UY=0. 

UZ=DIRZ/DIR 

C 

GO TO 2000 
C 

3000 IF(K.LT.l .QR*K.GT »KTOT »GR. 

+ RX*LT ♦ O.OR *RX.GT.XTOT)THEN 

PRINT** 'LEADER SEGMENT OUT OF BOUNDS' 

WRITE! 3f*> 'LEADER SEGMENT OUT OF BOUNDS' 

END IF 

IF < EF AMB . LT ♦ BDF ) THEN 

PRINT** 'LEADER SEGMENT ENTERED REGION LT BREAK DOWN' 
URITE<3»*) 'LEADER SEGMENT ENTERED REGION LT BREAK DOWN' 
END IF 
C 

CLOSE ( 1 ) 

CL0SEC3) 

STOP 

END 


SUBROUTINE RANDOM! I SEED* RAND) 
RANDOM NUMBER GENERATOR 

DOUBLE PRECISION E31fE32fSIfSJ 

E31=2***31 

E32=2.**32 

SI-DFLOAT ( ISEED) 

S J =MOD < 6 906? . *S 1 + 1 . r E32 ) 
IF(ABS(SJ> *LT .£31 ) THEN 
SI=SJ 

ELSEIF(SJ»GE*E31 )THEN 

Sl=SJ-E32 

ELSE 

SI=SJ+E32 
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END IF 

IF(SI *G£ .0. )THEN 

RAND*SI/E32 

ELSE 

RAND=1 .+SI/E32 
END IF 

ISEED*INT(SI ) 

RETURN 

END 


LDPM-INTERACTIVE 
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C INTERACTIVE VERSION OF LDPM 
C 

C PROGRAM DETERMINES DIRECTION OF NEXT LEADER SEGMENT DUE TO 
C ELECTRON WHICH WILL FORM AN IONIZED PATH TO 
C THE STEPPED LEADER SEGMENT TIP IN THE LEAST TIME 
C E FIELD TOTAL = E FIELD LEADER t E FIELD PREVIOUS LEADER 
C + AMBIENT E FIELD 

C 

COMMON/MAIN/FEX ( 9409 ) * FEZ ( 9409) 

COMMON/ VAR/NCX( 10) *NCZ( 10) *NCK( 10) *QN( 10 ) 

COMMON/DAT/NCELXf NCELZf ICX»ICZf NCfOP 
CHARACTER!!© INFILE *0FILE1 .0FILE3 
C 

QFILE1 =' LDPM .DAT ' ! BINARY PLOTFILE OF LIGHTNING PATH 

QFILE3=' LDPM. PRINT' ! PRINT VALUES OF LIGHTNING MODEL 
C 
C 

PRINT!* 'LIGHTNING PROPAGATION - TORTUOSITY MODEL' 

PRINT!* 'DEVELOPED BY E.M.A. - AUGUST 1986' 

PRINT!* ' 

PRINT!* 'SPECIFY GRID INPUT DATA' 

PRINT!* ' 

PRINT!* 'GRID SIZE - ENTER THE TOTAL NUMBER OF CELLS' 
PRINT!* 'EACH CELL IS 200* >: 200*' 

PRINT!* 'ALONG THE X (HORZ AXIS)* ALONG THE Z (VERT AXIS)' 
PRINT!* 'ENTER X*Z (MAX IS 97*97)' 

READ!* NCELX*NCELZ 
C 

PRINT!* ' 

PRINT!* 'ENTER THE INITIAL CELL COORD (X*Z) OF DISCHARGE' 
PRINT!* 'ENTER X*Z' 

READ!* ICX*ICZ 

PRINT!* 'ENTER THE CHARGE ( +.-COUL) ' 

READ!* OP 
C 

PRINT!. ' 

PRINT!* 'YOU MAY SPECIFY ANY CHARGE RE6I0N COORDINATE ' 
PRINT!* 'OTHER THAN THE INITIAL DISCHARGE COORDINATE* ' 
PRINT!. 'IF YOU DONT WANT THESE REGIONS ENTER 0 ' 

PRINT!* 'IT IS SUGGESTED THAT OVER THE ENTIRE GRID ' 
PRINT!* 'THE SUM OF THE TOTAL CHARGE * 0 ' 

PRINT!* 'ENTER THE ♦ OF CHARGE REGIONS (MAX IS 10)' 
PRINT!* '(NOT INCLUDING THE INITIAL DISCHARGE COORDINATED 
READ!* NC 

IF (NC.EQ *0)G0TQ 101 

PRINT!. 'ENTER THE X.Z PAIRS AND THtIR CHARGE(-f +COUL) ' 
PRINT!, 'ie: XI *Z1 .CHARGE1 ' 

PRINT!, ' X2.Z2.CHARGE2' 

DO 10 1=1 *NC 

10 READ!* NCX ( I ) * NCZ ( I ) * QN( I ) 

DO 11 1=1 *NC 

11 NCK ( I ) =NCX ( I ) +NCELX! ( NCZ ( I ) -1 ) 

C 

101 PRINT!* ' 

PRINT!* 'THE PROGRAM USES A MONTE CARLO METHOD TO' 

PRINT!* 'GENERATE FREE ELECTRONS AROUND EACH STEPPED' 
PRINT!* 'LEADER TIP. AN ASSUMED DENSITY IS *2elec/iit!3' 
PRINT!* '300 ELECTRONS IS RECOMMENDED' 

PRINT!, 'LESS THAN 100 WILL LIMIT TORTUOSITY' 

PRINT!* 'ENTER THE NUMBER OF ELECTRONS DESIRED' 


ORIGINAL PAGE IS 
OF POOR QUALITY 


ORIGINAL PAGE B 
BE POOR QUALTO5 


READ*; NE 
C 

PRINT*, ' 

PRINT*, 'THE MONTE CARLO ROUTINE REQUIRES A SEED INTEGER' 
PRINT#* 'ENTER ANY INTEGER UP TO 5 DIGITS I LT 32000*' 
READ#* I SEED 
C 

PRINT*, ' 

PRINT*, 'PROGRAM GENERATES A PRINT FILE OF THE MODEL' 
PRINT*, 'VALUES S A BINARY PLOTFILE IN AN XrZ OUTPUT' 
PRINT*, 'AMBIENT E FIELDS ARE PRINTED OUT' 

PRINT*, 'THESE FILES ARE LDPN.PRINT, LDPM.DAT, AMB.DAT' 
PRINT*, ' 

PRINT*, tttttXttttttttttttttttttttttttttttttttt' 

PRINT*, ' 

C 

C 

C CALCULATE AMBIENT E FIELDS 
CALL EGRID 

C INPUT VARIABLES 

C CL S LEADER LENGTH ELAMBDA=LEADER CHARGE DENSITY 

C BDF=BREAK DOWN FIELD EK=l/4PIE*so NL-LEADER t 

C DZ=Z CELL HEIGHT DX=X CELL WIDTH 

C NE=* FREE ELECTRONS ED=ELECTRON DENSITY <.2/»**3) 

C Q=MAX ALLOWABLE ANGLE BTU EFIELD l DIRECTION VECTORS 
C 

DATA CL , ELAMBDA , BDF/50 .,.001,1 .5E5/ 

DATA PI , EK »NL/3.1 415927, 9*E9,0/ 

DATA DZ,DX,NE, ED, Q/200.,200.,500,. 2, 12. 5/ 

DATA EXM1 ,EYMl ,EZMl/3*0./ 

CONST^EKtELAMBDA 
VOL=FLOAT(NE)/EIi 
RADIUS^ ( .75*V0L/PI >*<< 1 ./3. ) 

KTOT=NCELX*NCELZ 

XTOT=NCELX*DX 

C 

C INITIALIZE LEADER PROPAGATION 
C INPUT DATA: Rl, AMBIENT E FIELD, U1 
C COORDINATES OF INITIAL POINT ON GRID PROBLEM SPACE 
RX=DX*(FL0AT(ICX)-.5> 

RY-O. 

RZ= r DZ*(FL0AT(ICZ)-.5) 

C 

IC*ICX+NC£LX*(ICZ-1) 

ET0T=SQRT(FEX(IC)**2+FEZ(IC)**2) 

C DEFINE UNIT VECTOR 
UX=FEX( IO/ETQT 
UY-O. 

UZ=FEZ( IO/ETQT 
C 

0PEN(1 ,FILE-0FILE1, FORM* ' UNFORMATTED ' , STATUS* ' NEW ' ) 
GPEN<3,FIL£-0FILE3, STATUS*'N£W' ) 

C 

C START LOOP 
2000 NL*NL+1 
C 

C DETERMINE CELL POSITION OF NEW LEADER SEGMENT TIP ON OVERLAY 

C USE APPROPRIATE AMBIENT E FIELD MAGNITUDE AND DIRECTION 

C K IS THE CELL NUMBER OF THE LEADER SEGMENT TIP LOCATION 
KX=ANINT < (RX/DX) + .5) 
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KZ=ANINT ( (RZ/DZH»5) 

K=KXfNCELX*(KZ-l) 

EFAMB=SQRT(FEX(K) M2+FEZ ( K ) **2 ) 

C 

C TERMINATE IF 

C LEADER HAS PROPAGATED OUTSIDE THE GRID 
C LEADER ENTERING CELL WITH AN AMBIENT FIELD LT 1.5M V/M 
C LEADER ENTERS INTO CHARGED REGION 
IF(K*LT.l*OR*K.GT *KTQT -OR* 

+ RX.LT .O.OR *RX *GT *XTOT )GOTO 3000 
C 

URITE(3>*) 'LEADER SEGMENT NUMBER 'iNL 

WRITE(3# *) 'CELL NUMBER '#K#' <X#Z> COORD'#KX#KZ 

URITE<3#*)'T0T E FIELD V/b'#E TOT#' AMB E FIELD '#EFAMB 

URITE(3# *) ' AMB EX 'rFEXOO#' AMB EZ 'fFEZ(K) 

URITE(3»*) ' X iet» froi 0#0'#RX#' Z Bet. froB 0#0 '#RZ 

URITE<3#*>' 

PRINT** 'LEADER SEGMENT NUMBER ' »NL 
PRINT*! 'CELL NUMBER '#K#' <X#Z> COORD'#KX»KZ 

PRINT*# 'TOT E FIELD V/b'fETQT#' AMB E FIELD ' fEFAMB 
PRINT*# 'AMB EX SFEX(K)#' AMB EZ 'fFEZ(K) 

PRINT*# 'X aet* froB 0#0'»RX#' Z Bet. froa 0#0 '#RZ 
PRINT*# ' 

URITE( 1 )RX#RZ 
C 

IF(EFAMB*LT.BDF)GOTO 3000 
IF(NC.EQ.O.)GOTO 55 
DO 111 1 1 = 1 » NC 

111 IF( K »EQ*NCK( 1 1 ) )GOTO 3000 

SHORTEST TIME OF ELECTRON TO FORM IONIZATION CHANNEL 
C THIS TIME DETERMINES DIRECTION OF NEW LEADER SEGMENT 
55 TMIN=1*E20 

C 

C USE MONTE CARLO METHOD TO GENERATE RANDOM FREE ELECTRONS 
C CALCULATE ELECTRON POSITIONS NEAR LEADER TIP 
DO 100 1=1 #NE 

C 

C IN WORLD COORDINATES 

CALL RANDOM< ISEEDfRAND) 

RE=RAND*RADIUS 

C RE=SQRT(RAND*(RADIUS*<2) ) 

CALL RANDOM(ISEED#RAND) 

th= ( rand ) *pi 

C TH=AC0S((RAND*2)-1) 

CALL RANDOH< ISEEDfRAND) 

PH=(RAND)*2.*PI 

X=RE IS IN ( TH > *COS ( PH ) +RX+UX*CL 

Y=RE *SIN< TH ) *SIN( PH ) +RY+UY*CL 

Z=RE*COS(TH)FRZ+UZ*CL 

F!HX=X-RX-UX*CL 

RMY=Y-RY-UY*CL 

RMZ=Z-RZ-UZ*CL 

RM=SORT ( RMX**2+RMY**2+RMZ**2 ) 

C 

C COMPUTE E FIELD COMPONENTS FOR Ith ELECTRON 
AX=X-RX 
AY=Y-RY 
AZ=Z-RZ 

A=SORT ( AX**2+AY**2+AZ**2 ) 
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IF < A * LT « 1 *E-6)G0TQ 100 

B=AX*UX+AY*UY+AZ*UZ 

AS=A**2 

BS=B**2 

IF ( AS-BS*LT . 1 *E-6)G0T0 100 
P=CONST/(AS-0$> 

D=CL**2-2>tBICL+AS 
IF(D*LE.O. )G0T0 100 
D=S0RT'D) 

C 

C E FIELD COMPONENTS DUE TO LEADER SEGMENT AND AMBIENT FIELD 
EX=P*< ( (AX*(CL-B)+UX* <AS-B*CL> )/D)- 
+ { (UX*AS-AX*B)/A) )+FEX(K) 

EY*Pt< ( ( AY*(CL-B)+UY*<AS-B*CL) >/D>- 
+ ( (UY*AS-AY#B)/A) ) l+FEYOO 

EZ=P<( < < AZ*(CL-BHUZ*(AS-B*CL> >/D>- 
+ ( (UZ*AS-AZ*B)/A) )+FEZ(K> 

C 

IF(NL.LE.1)G0T0 1000 
AXM1=X-RXM1 
AYM1=Y-RYM1 
AZM1=Z-RZM1 

AMI =SQRT ( AXMl*«2+AYMU«2+AZMltt2 ) 

IF (AMI *LT* 1 .E~A)G0TG 100 

BM1=AXM1*UXMHAYH1*UYMHAZM1*UZM1 

ASM1=AM1**2 

BSM1-BM1I*2 

IF(ASM1-BSM1 *LT. 1 .E-6) GOTO 100 
PM1=CQNST/(ASM1-BSM1 ) 

DH1=CL«2-2.*BM1*CL+ASH1 
IF(DM1 .LE*0. JGQTO 100 
DM1=SQRT ( DM1 ) 

C 

C E FIELD COMPONENTS DUE TO PREVIOUS LEADER SEGMENT 

EXM1=PM1*( ( ( AXM1*(CL-BM1 )+UXMlt(ASMl-BMl*CL) )/DMl>- 
+ ( (UXMlf ASMl'AXMKBMl )/AMl ) ) 

EYMl*PMl»(((AYNU(CL-BMmUYMl*(ASMl-BMUCL))/DMl)- 
+ ( (UYN1*ASM1-AYM1<BM1)/AH1) ) 

EZM1-PM1*( ( (AZMI KCL-BM1 >+UZHlt(ASMl-BMltCL) )/DMl )- 
+ ( ( UZMltASMl -AZMI *BM1 ) /AMI) ) 

C 

C TOTAL E FIELD COMPONENTS FOR Ith ELECTRON 
1000 EX=EX+EXM1 
EY=EY+EYM1 
EZ=EZFEZ«1 

ET=saRT(EX**2+EY**2+EZt*2>+.01 

DL=(RMX*EX+RMY*EYFRMZtEZ)/(RMfET) 

C 

C CHOOSE APPROPRIATE ELECTRON TO FORM IONIZED CHANNEL 
C TO THE LEADER TIP IN THE SHORTEST TIME 

IF<ET.LT.<3.0E6)*AND.(RH/ET).LT.THIN> 

+ AND*ACQS<ABS(DL)).LE.Q*PI/180.)THEN 
TMIN=RM/ET 
ETOT=ET 
DIRX=EX 
DIRY=E < 

DIRZ-EZ 
END IF 

ICO CONTINUE 
C 
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C CALCULATE DIRECTION FOR NEXT LEADER SEGMENT 
DRX=CL*UX 
DRY=CL!UY 
DRZ=CL*UZ 
RXM1=RX 
RYM1=RY 
RZMZ=RZ 
RX=RX+DRX 
RY=RY+DRY 
RZ-RZFDRZ 
C 

C NEU UNIT VECTOR 

C DIR*SQRT(DIRX!t2+DIRY**2+DIRZtt2) 

DIR*SGRT(DIRXtt2+DIRZW2> ! PRQJ'N ONTO 2D GRID 
IF(DIR*LT»1 *E-6)GOTO 2000 
UXM1=UX 
UYM1=UY 
UZM1=UZ 
UX=DIRX/DIR 
C UY=DIRY/DIR 

UY=0* 

UZ=DIRZ/DIR 

C 

GO TO 2000 
C 

3000 IF(K*LT. 1 *OR*K.GT #KTOT*OR* 

+ RX.LT.O*OR*RX.GT*XTOT) THEN 
PRINT!* 'LEADER OUT OF BOUNDS' 

WRITE<3*!> 'LEADER OUT OF BOUNDS ' 

END IF 

IF(EFAMB.LT*BDF)THEN 

PRINT! » 'LEADER ENTERED REGION LT BREAK DOWN' 
WRITE (3* *> 'LEADER ENTERED REGION LT BREAK DOWN' 
END IF 

DO 112 II 3 1 >NC 

IF(K.EQ.NCKGI) )THEN 

PRINTIf 'LEADER ENTERED CHARGED REGION' 

URITE(3* t) 'LEADER ENTERED CHARGED REGION' 

END IF 

112 CONTINUE 
C 

CLOSE ( 1 ) 

CLGSE<3> 

STOP 

END 

C 

C 

SUBROUTINE RANDOM ( ISEEDf RAND) 

C RANDOM NUMBER GENERATOR 
C 

DOUBLE PRECISION E31 *E32?SI »SJ 

E31=2*t!31 

E32=2.!!32 

SI=DFLOAT ( ISEED ) 

SJ=MQD< 69069 . tSI+1 • fE32) 

IF(ABS(SJKLT>E31)THEN 

SI=SJ 

ELSEIF(SJ.GE.E31 )THEN 

SI=SJ-£32 

ELSE 
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SI=SJ+E32 
END IF 

IF(SI *GE .0. )THEN 

RAND*SI/E32 

ELSE 

RAND=1 • FSI/E32 
END IF 

ISEED=INT(SI) 

RETURN 

END 

C 

C 

SUBROUTINE EGRID 


original 

0£ POOR 
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C CALCULATION OF AMBIENT ELECTRIC FIELDS 
C 

DIMENSION X(?409)»Z<9409) 

DIMENSION XMRN<5)t B(5) 

COMMON/ MA IN/FEX ( 9409 ) , FEZ (9409 ) 

COMMON/ VAR/NCX( 10) rNCZ(10) rNCK( 10) »QN( 10) 
COMMON/DAT /NCELX >NCELZ » ICX» ICZ»NC»GP 
DATA DX » DZ r DP *P 1/200* *200* t9*E9r 3. 1415927/ 
C 

OPEN ( 4 , F I LE= ' AMB * DAT ' i STATUS* 'NEW ' ) 

C 

DO 10 J=1 »NCELZ 
DO 10 1=1 f NCELX 
K*I+NCELXf < J-l ) 

X(K)=DXt( FLOAT ( I)-*5) 

Z(K)=DZt (FLOAT( J)-*5) 

K1=K 

10 CONTINUE 


C 

LP*ICX+NCELX*(ICZ-1> 

C CALCULATE RESULTANT E FIELDS FROM CHARGE CENTERS 
DO 20 K=1*K1 
IF(NC*EQ.O)THEN 

IF(X(K) *EQ. X(LP) *AND*Z(K) *EQ*Z(LP> )GOTO 20 
END IF 

IF(NC*E0*0)G0T0 23 
DO 333 1=1 ,NC 

IF < X < K ) *EQ.X(LP) .AND.Z(K) *E0»Z(LP ) * OR ♦ 

+ X(K) *EQ*X(NCK(I)) * AND. Z(K) *EQ* Z(NCK ( I ) ) )GOTO 20 
333 CONTINUE 

C 

23 XMRP=SQRT( (X(K)-X(LP) )M2+(Z(K>-Z(LP) >W2) 
A*DP*QP/XMRP«3 

FEX(K )*0. 

FEZ(K)=0. 

IF(NC*E0.0)G0T0 24 
DO 444 1=1 ? NC 

XMRN< I )=SGRT ( <X(K)-X(NCK(I>) )$*2+<Z(K)-Z<NCK(I > >><*2> 
B ( I ) =DP*QN < I ) /XMRN ( I ) ttZ 
FEX(K)=FEX(K)+B(I)*(X(K)-X(NCK(I>)> 
FEZ(K)=FEZ(K)+B(I)XZ(K)-Z(NCK(I))> 

444 CONTINUE 

C 

24 FEX(K)=FEX(K)4A*(X(K)-X(LP)> 
FEZ(K)=FEZ<KHA*<Z(K>-Z<LP>> 

20 CONTINUE 

FEX(LP)*FEX(LP+1 ) 
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FEZ(IP)=FEZCLPM> 

IF(NC.EQ*0>GOTO 25 

DO 555 1*1 >NC 

FEX(NCK< I ) )=FEX(NCK( I)-l ) 

FEZ(NCK( I) )*FEZ<NCK< I)-l) 

555 CONTINUE 

25 DO 26 K*1fK1 

URITE(4t f )K»FEX(K) »FEZ(K)»X(K) t Z(K) 

26 CONTINUE 
CLQSE<4> 

RETURN 

END 
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1.0 


INTRODUCTION 


The maximum observed electric fields in thunderstorms are in general 
smaller than is necessary to initiate corona or sparking that can lead to lightning 
events. Field magnitudes in the range of 1 to 2 megavolts per meter are necessary in 
pure air at typical thunderstorm altitudes, whereas maximum measured fields are of 
the order of 400 kilovolts per meter. Hence it is necessary to look for other 
mechanisms which can locally raise the field strength to the necessary intensity. 
Thunderstorm particles are one source of local field enhancement which may provide 
the necessary mechanism. These particles (cloud droplets, rain, ice, hail, snow, 
graupel), although not having large conductivities, do have small electrical relaxation 
times on the scale of non-lightning electric field changes in a cloud. This implies that 
the particles act as perfect conductors in a thundercloud, and are able to locally 
enhance the ambient electric field. They are also able to hold electrical charges which 
can alter the way in which field enhancement occurs. 

The purpose of the research presented here is to attempt to determine what 
effect thunderstorm particles have on the lightning initiation process. Are there 
preferred areas within the storm where particular types of particles exist that make 
lightning initiation more likely? Are large particles or large particle densities more 
effective at producing electrical breakdown in a particular volume? To answer these 
questions a model has been developed to roughtly calculate the ambient field 
necessary to cause air breakdown in a thunderstorm in the presence of ensembles of 
particles. The size, shape, and density of the particles are variable inputs to the 
model, as is the altitude at which the model does its calculation. The altitude (or air 
density) is the most significant parameter in calculating the pure air electron avalanche 
rate [1], and the particle parameters determine how that pure air rate is altered. Using 
this model a table can be constructed giving the nominal air breakdown field as a 
function of the above parameters. The table can then be used in conjunction with 
large scale thunderstorm models [2] which calculate particle types and densities as a 
function of location in a cloud to predict the most likely regions for lightning initiation. 

Electrical breakdown of air occurs because of a growth in the electron 
density in the air at a given location. This occurs when electrons are accelerated in an 
ambient electric field to energies large enough to ionize air molecules through 
collisions. The process is known as an electron "avalanche." The growth continues 
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approximately exponentially until the electron density is large enough to produce 
currents which neutralize the ambient field locally. These electron currents, while 
neutralizing a local electric field, also raise the field at other nearby locations, possibly 
causing the electron avalanche to propagate away from its initial location. In this way 
a lightning channel may be initiated. 

This report investigates specifically the way in which the electron avalanche 
rate varies as a function of particle parameters, air density, and electric field. Nominal 
air breakdown fields are predicted for these variables by calculating the smallest 
ambient field which will produce sufficient electron density to lower that field. This 
involves the balancing of the electron avalanche rate with the primary loss mechanism 
for electrons, attachment to neutrals to form negative ions. Hence the nominal air 
breakdown field is that for which the avalanche rate equals the electron attachment 
rate. 


The results, to be presented in greater detail later, show that large particles, 
such as snowflakes, are the most effective at altering the nominal air breakdown field, 
and volumes containing these particles are the most likely locations for lightning 
initiation. Small particles, even if large densities are present, are not nearly as 
effective in altering the nominal air breakdown field. 

2.0 CORONA GROWTH AROUND PARTICLES IN A UNIFORM 

ELECTRIC FIELD 

The work presented in this section has been partially published in a 
previous report [3] and is also presented here in the interest of completeness. 

By locally redistributing the electric field energy density in a region of space, 
thunderstorm particles can affect the way in which electrons and ions are produced in 
an electron avalanche. In particular, a uniform electric field which is below breakdown 
intensity may be locally enhanced to a level which produces corona around the 
particles. This occurs because the avalanche rate is strongly nonlinear with electric 
field, so that small electric field changes can result in large changes in charged 
particle production. 
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As an example of this consider the case of a spherical particle in a uniform 
static electric field, E 0) with magnitude ninety percent of air breakdown intensity. This 
is shown schematically in Figure 2.1. The electrical conductivity of the particle is 
important only for very fast field changes (of the order of microseconds), and in general 
the particle acts as a perfect conductor on the time scale of large scale thunderstorm 
fields. Because of this the well-known analytic development for the altered field 
distribution around a conducting sphere in a uniform electric field may be used. Figure 
2.1 shows the regions around the particle in which fields larger and smaller than the 
ambient field exist. In addition the regions in which the field magnitude exceeds air 
breakdown level, Et>, are shown. These latter regions depend on the magnitude of the 
ambient field with respect to breakdown field. That is, if the ambient field drops, the 
breakdown region will shrink, finally disappearing altogether when the ambient field 
goes below one-third of the breakdown field. The other regions, those of enhanced 
and decreased field, are not dependent on ambient field strength. 

The preceding discussion was based on a spherical particle. The situation 
for more irregular particles such as ice crystals or snow flakes is more complicated. 
Analytic solutions for these shapes do not exist, so approximations must be made in 
specifying the local field distribution around the particle. In general, however, a 
polarized particle will behave like an electric dipole to some level of approximation. 
The assumption used in the development presented here is to treat the local field 
distribution as the same as that around a spherical particle, but to raise the maximum 
field enhancement to that determined from a previous study [4] (e.g., enhancement of 
five for columnar ice crystals and nine for a typical snowflake). 

The importance of Figure 2.1 is that it shows that a conducting particle can 
cause electron avalanching to occur in a region of space which would have no 
avalanching in the particle’s absence. The details of how much avalanching occurs 
because of the particle and the effects of irregularly shaped particles and particle 
density are left for the mathematical development of the next section. 
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2.1 


Mathematical Development 


To facilitate the mathematical development to follow, some symbols will first 
be defined. The development will be done in a spherical coordinate system, assuming 
no azimuthal variation. The ambient field is assumed to be oriented along the polar 
axis. The assumption of no azimuthal variation is strictly true only for spherical 
particles or cylindrical particles oriented with their axis along the polar axis. Particles 
not satisfying this criterion will be included in an approximate sense. 

Table 2.1 

Definition of Symbols Used in Section 2.3 


a 

- approximate half size of particle (radius in the 
case of a spherical particle) 

e 0 

- ambient electric field intensity 

E b 

- breakdown electric field intensity 

E r 

- radial electric field intensity 

E e 

- theta component of electric field intensity 

E z 

- component of electric field intensity along polar 
axis 

^max 

- maximum electric field enhancement factor 
(particle dependent) 

n 

- density of thunderstorm particles 

N 

- number of thunderstorm particles in a finite 
difference cell 

P 

- air density relative to sea level density 

AV 

- volume around a particle in which field is above E b 




The objective of the analysis to follow is to develop a modified electron 
avalanche rate which accounts for the presence of thunderstorm particles. Previous 
work has assumed that corona formation occurs in clear air only. To include particles, 
it is necessary first to calculate the avalanche occurring around a single particle and 
then to consider the effects of a particle density in a given region of space. The electric 
field around a single particle is investigated first. 

The field distribution outside of a spherical conducting particle placed in a 
uniform electric field is given by, 

E r (r,9) = E 0 (l + 2 a jj)cose (2.1) 

E 6 ( r .6) = -E 0 (l - 4j- ) sin 6 

Figure 2.2 shows the spherical coordinate system used in the analysis. For use here 
this expression is generalized to approximate the field distribution around irregularly 
shaped particles as follows, 


E r (r.e> = Ejl +(f max -1)-^ ] cose (2.2) 

E e (r.e> = -E 0 (l--*r )*"« 

Note that the only difference between Equations 2.1 and 2.2)is that the 2 in the first of 
Equations 2.1 has been replaced by the term f max -1 in Equation 2.2. This reduces to 

Equation 2.1 in the case of a sphere (for which f max = 3). The justification for using 

Equation 2.2 in the present model is not rigorous. Any conducting particle placed in a 
uniform electric field must polarize, and the resulting field distribution must 
in some sense resemble a dipole field. The assumption is made here, mainly for the 
purpose of making further analytic progress on particles in general, that the distribution 
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Figure 2.2 


Spherical Coordinate System Used in Analysis of 
Thunderstorm Particles 
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is similar to that around a spherical particle, but with differing maximum field 
enhancements. 

In addition, the electric field used in the calculations to follow is the 
component along the polar axis of the coordinate system. Strictly speaking, the total 
electric field is more appropriate. However, the form of the polar component is much 
easier to work with analytically, and introduces only a small error. The error is small 
because the axial oriented field is a good approximation to the total field at most 
locations. The axial field is easily calculated from Equation 2.2 and is shown below. 

E z (r.9) = E 0 [l + (f max cos 2 e -D-J-] (2.3) 


This expression involves three physical quantities of interest, the ambient 
electric field, the size of the particle, and its shape (included in f max ). These three 

quantities, along with the breakdown field strength at a particular altitude and the 
particle density, completely determine the altered avalanche rate from the clear air 

expression. The next step in the analysis is to specify the volume around the particle 
over which avalanching occurs. In other words, over what volume is E z larger than 

E b ? The condition for this is expressed in equation (2.4). 

< 2 - 4 > 


There are three possible regimes for Equation 2.4. If E Q is larger than E b , 

then Equation 2.4 is satisfied in all of space. This is of course the trivial case in which 
air breakdown would occur in the absence of particles. A second regime is one in 
which Equation 2.4 is satisfied nowhere in space. This occurs when f max is not large 
enough to raise the maximum enhanced field above E fa . In this case air breakdown 

is absent even in the presence of particles. The middle regime in which Equation 2.4 

is satisfied in a limited region of space is the most interesting. The conditions 

fmax ^ j^and < define this regime. Here air breakdown occurs when particles 
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are present but not when they are absent. From Equation 2.4 the volume over which 
this occurs can be specified as below. 


0 < 0 < cos 




max^o 


m A 




ax 


COS 


W^T 




= Ba 


(2.5) 


Equation 2.5 represents a two dimensional projection of a three dimensional 
volume. The full volume is obtained by rotating the area specified in Equation 2.5 
around the polar axis of the particle. Note also that there are two of these volumes, 
one on either end of a symmetrically shaped particle. An expression for the volume 
over which the enhanced electric field is above E can now be derived by integrating 
over the limits of Equation 2.5 and multiplying by two to account for both ends of the 
particle. Formally this is written, 



o o a 


(2.6) 


The details of the integration are straightforward, and the final expression for 
the avalanche volume is, 


AV = 


4rca v- 


^o^max 

E. - E. 


1 +2 


E o *max 


\3/2 

■J - 3 


E o Vnax 


■)] 


(2.7) 


The most important thing to notice about Equation 2.7 is the way in which AV 
depends on various quantities. A large AV implies that the particle has a significant 
effect on electron avalanching, and a small AV implies very little effect. Therefore it is 
clear that larger particles (large a) are most effective at altering the avalanche rate, 
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since AV varies as the cube of the particle dimension. The shape is also involved, but 
is less important, in the sense that AV is nearly linearly dependent on f__ v . Of course it 
must be remembered that f max is quite important in determining whether there will be 

any avalanching at all around the particle. 

Equation 2.7 gives the volume around a single particle in which avalanching 
will occur. To implement this concept in a numerical model, it is necessary to include 
the density of particles which is present in the region of interest. Hence it is necessary 
to calculate how much of a given volume has fields above breakdown level. The total 
volume with field above breakdown strength in a given volume is then, 

AV t = NAV (2.8) 

where N is the number of particles in the given volume and is numerically equal to the 
particle number density multiplied by that volume. In Equation 2.8 the assumption has 
been made that the particles are noninteracting in the sense that none of the individual 
enhanced field volumes overlap. This appears to be a reasonable approximation for 
typical particle densities in thunderstorms. 

The next step in the analysis is to calculate the number of electron-ion pairs 
produced in AV^ per unit time. This is essentially the desired avalanche rate from 

particles. The rigorous method to find this quantity is to integrate the (strongly electric 
field dependent) avalanche rate over AV. This is computationally inefficient for large 
scale air breakdown codes, however, so the method used here is somewhat different. 
The alternate method is to define an average electric field in the breakdown volume 
and to simply evaluate the altered avalanche rate for that electric field. One possible 
choice of this average field is to set it equal to the arithmetic average of the smallest 
and largest fields in the volume. These fields are E b and f max E 0 , respectively. 

However, because the largest field exists only in a very small subvolume of AV it 
would be better to favor the lower fields in any average. Therefore the geometric 
mean is used, as shown below. 


Average field = VE„fmax E o 


(2.9) 
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The geometric mean is always less than or equal to the arithmetic mean, so this 
should be a more appropriate value to use in the model. 

Finally the avalanche rate from particles may be formally written, 

AV 

^particles = G pure air( avera 9 e ,ield >-V~ < 2 - 10 > 

The last matter to be considered is the issue of irregularly shaped particles 
which may not be aligned with their axis along the ambient electric field. Although 
electrostatic forces would tend to align these particles, it is likely that hydrodynamic 
forces from winds and small scale turbulence would overwhelm the electrostatic force. 
Hence the particles are probably randomly oriented with respect to the field direction. 
This is taken into account in the present model by assuming that on the average half of 
the particles are aligned along the field and half are aligned across the field. The 
particles across the field make no contribution to the particle dependent avalanche 
rate. In the model this is accomplished simply by dividing the true particle density by 

two to arrive at an effective density. It should be understood that this is done only for 
the case of nonspherical particles (f max ^3), because for spherical particles alignment 

is unimportant. 

3.0 RESULTS 

The immediate result which can be calculated from the thunderstorm particle 
model is the effective ambient field which causes air breakdown to occur. This should 
be lower in the presence of thunderstorm particles than in their absence. An 
approximation to the breakdown field is achieved by finding the field level at which the 
avalanche rate becomes larger than the electron attachment rate. Although this is not 
a strict criterion for air breakdown it does serve to illustrate the difference in corona 
growth caused by particles. Table 3.1 shows the nominal air breakdown field at 
various altitudes for a number of typical thunderstorm particle environments. 

It should be noted that not all altitudes shown in Table 3.1 are appropriate 
for each of the environments. For instance, the ice particles and snowflakes of Cases 
6-14 would not be expected below the freezing level in a thunderstorm. 
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Table 3.1 

Nominal Air Breakdown Field (MV/m) as a Function of Air Density and 

Particle Parameters 

(Note: Not All Air Densities Are Appropriate for Particular Particles) 


\ p 
Castf\^ 

.3 

.4 

.5 

.6 

.7 

.8 

.9 

1.0 

1 

.64 

.86 

1.07 

1.29 

1.50 

1.72 

1.93 

2.15 

2 

.64 

.86 

1.07 

1.29 

1.50 

1.72 

1.93 

2.15 

3 

.64 

.86 

1.07 

1.29 

1.50 

1.72 

1.93 

2.15 

4 

.64 

.86 

1.07 

1.29 

1.50 

1.72 

1.93 

2.15 

5 

.64 

.86 

1.07 

1.29 

1.50 

1.72 

1.93 

2.15 

6 

.64 

.86 

1.07 

1.29 

1.50 

1.72 

1.93 

2.15 

7 

.64 

.96 

1.07 

1.29 

1.50 

1.72 

1.93 

2.15 

8 

.64 

.86 

1.07 

1.29 

1.50 

1.72 

1.93 

2.15 

9 

.64 

.86 

1.07 

1.29 

1.50 

1.72 

1.93 

2.15 

10 

.64 

.86 

1.07 

1.29 

1.50 

1.72 

1.93 

2.15 

1 1 

.64 

.86 

1.07 

1.29 

1.50 

1.72 

1.93 

2.15 

12 

.64 

.86 

1.07 

1.29 

1.50 

1.72 

1.93 

2.15 

13 

.64 

.85 

1.07 

1.28 

1.50 

1.71 

1.93 

2.14 

14 

.62 

.83 

1.04 

1.25 

1.45 

1.66 

1.87 

2.08 


p = relative air density 
Key to cases: 


1 - 

pure air 

8 

2 - 

r = 2.5 p, n = 10 8 nr 3 , f max = 3 

9 

3 - 

r = 10 p, n = 10 8 nv 3 , f max = 3 

10 

4 - 

r = 2.5 p, n = 10 9 nr 3 , fmax = 3 

11 

5 - 

r=10p.n = 109m-3,f max = 3 

12 

6- 

r = 50 p, n = 10 4 m' 3 , f max = 5 

13 

7 - 

r = 500 p, n =1 0 4 nr 3 , f max = 5 

14 


r = 50 p, n = 1 0 5 nr 3 , f max = 9 
r = 500 p, n = 10 5 nr 3 , fmax = 9 


r = particle radius 
n = particle density 
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Although the contents of Table 3.1 are rather dull, they illustrate the principle 
that small particles do not significantly change the air breakdown field, even when 
present at high densities. Only the very largest particles with large enhancement 
factors, such as snowflakes, are able to change the breakdown field by more than 1%. 
Hence, if e ectric fields everywhere in the cloud were equal, it is slightly more likely 
that regions having this type of particle will initiate lightning. However, the field 
variation throughout a typical cloud will certainly swamp the small selection effect 
toward large particles shown in Table 3.1, so in practice the large particle regions are 
probably no more or less likely to initiate lightning than any other regions. 

4.0 CHARGED PARTICLE CONSIDERATIONS 

The previous analysis assumed that the thunderstorm particles were 
uncharged. A net charge on the particles introduces several complicating factors. The 
total electric field around the particle can be represented as the sum of the field from 
the analysis of Section 2.1 and the field from the net charge. For spherical particles 
the field from the charge is just the well known field from a point charge. For more 
complicated shapes such as ice columns and snowflakes the field distribution 
becomes much more complex. For instance, for an ice column with a net charge the 
field is likely to be more characteristic of two point charges separated by the length of 
the column. For snowflakes sharp points are likely to behave as point charges, each 
comprising a fraction of the total charge. Hence the field from a net charge around a 
point where corona may occur may be roughly characterized as that around a point 
charge carrying a fraction of the total charge on the particle. That fraction may be 
estimated by the number of potential corona points on the particle: one for spheres, 
two for columns, and possibly six for undamaged snowflakes. 

The maximum total charge on the object is also a function of the type of 
particle and its size. Corona growth bleeds charge into the surrounding atmosphere 
when the total charge is large enough to produce air breakdown level fields around 
the particle. Clearly larger particles can hold larger charges, and spheres, because of 
their low enhancement factor, can hold larger charges than equally sized particles of 
other types. The upper limit on total charge for a particle is then the quantity which 
produces breakdown around a sphere of roughly equal size, as given below. 

Q < 47t£oa 2 Eb (4.1) 
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Here Eb is the breakdown field at the altitude of interest, and a the nominal radius of 
the particle. This is roughly 2 picocoulombs for a particle of radius 100 |i at an altitude 
for which the air breakdown field is 1 .5 megavolts per meter. 

The interaction of a net charge and an ambient field on the possible air 
breakdown around the particle is interesting. An uncharged particle has volumes of 
equal electric field levels on either side as was shown in Figure 2.1. When a net 
charge of either sign is placed on the particle, the field on one side of the particle is 
lowered and on the other side elevated. Qualitatively, this tends to suppress potential 
air breakdown on one side and enhance it on the other. While this effect reduces the 
net volume over which corona occurs by roughly a factor of two, the ambient field 
required to initiate corona on the enhanced side can be significantly reduced. This 
can be seen by writing down the total field component in the axial direction, in analogy 
with Equation 2.3. 

Ez(r,6) - E 0 [ 1 + (*max COS20 -1 )£ ] + COS8 (4.2) 

Here qf is the absolute value of the fraction of the total charge on the particle which 
contributes to the field<at the point in question. Unfortunately this expression is no 
longer a ve7 good approximation to the total field in the enhanced volume as in 
Equation 2.3, because of the spherical nature of the field from the charge. Because of 
this and the relatively complicated nature of Equation 4.2, a detailed analysis of the 
breakdown volume as was done in Section 2.1 is beyond the scope of this effort. 
However, a qualitative feel for the reduction in ambient field necessary to cause 
breakdown can be gained by investigating the maximum fields in Equation 4.2 in the 
cases of a particle with and without a net charge. The maximum field in each case 
occurs at 0 = 0, and r = a. Equation 4.2 becomes for this case, 


E"f‘ (a.o) = Eofmax + (4.3) 

Next assume that this maximum field is exactly the breakdown field Et>. The 
ambient field needed to produce this breakdown field is given by 


C - 1 5 
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E 0 = 
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Eb - 


q< 

47ie 0 a 2 


Then the ratio of the necessary ambient field needed to produce a breakdown in the 
charged and uncharged case is given by, 


E 0 charged/E 0 uncharged _ -j . 


a 

47ie 0 a 2 Eb 


(4.5) 


Hence, as one would expect, the reduction in ambient field necessary to produce air 
breakdown is linearly dependent on the net charge of the particle. 

This analysis, admittedly crude, can be used in conjunction with Table 3.1 to 
estimate the ambient field necessary to initiate lightning in a region of particles having 
net charge. 

5.0 EXPERIMENTAL INITIATION CRITERION 

In addition to the requirement that the absolute electric field magnitude 
exceed air breakdown in a volume for lightning initiation to occur, there is some 
experimental evidence, documented by Ruhnke [5], of a ceiling on the allowable field 
gradient. For field gradients above an absolute magnitude of 2.5 x 10 8 V/m 2 only 
corona will occur and no spark or lightning leader will be initiated. Field gradients 
below this level allow the initiation of a spark which may grow into a full leader. It is 
therefore of interest to evaluate the expressions for the electric field distribution around 
thunderstorm particles to determine whether this criterion is satisfied. 

The expression for the axially oriented field around a particle was derived in 
Equation 4.2. 


E 2 (r,0) 


0 

Eq [ 1 + (fmax cos 2 0 -1 )f] 


+ — sl_ 

47re 0 r 2 


COS0 


(5.1) 
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The largest electron avalanche and therefore the largest corona growth will occur at 
the point where the largest electric field exists. This is at the 0 = 0, r = a point. This is 
also the point at which the largest electric field gradient exists. Hence to satisfy the 
above criterion in the entire breakdown volume, it must be satisifed at this point. 
Because propagation of any potential spark will occur in the direction of the field, it is 
the radial derivative which is of most interest. This derivative is shown below. 


9Ez 

dr 


-T 2 (WCOS20-1) 


£ 

r 3 


9 L_ 

2jte 0 r 3 


COS0 


(5.2) 


Evaluating this expression at the point of maximum field gives, 


IT (r=a> 0=O) = ( f max * 1 ) * 2jt ^ a3 (5 ‘ 3) 

Note here that E 0 , qt, and f m ax all contribute in a linear fashion to the 
gradient. Hence particles with low net charges and small enhancement factors would 
appear to be more likely to satisfy the criterion. However, this is not really the case, 
because a decrease in qi on a particle results in a corresponding increase in E 0 in 
order to achieve a field of breakdown magnitude. The same is true of f ma x- E 0 ' s n0 * 
an independent variable in the equation. In effect Equation 5.3 can be evaluated with 
qt = 0, and f m ax = 3 to set order of magnitude values for field gradients of typical 
thunderstorm particles. 

The size of the particle, a, enters Equation 5.3 in an inverse fashion. Hence 
larger particles are more likely to satisfy the maximum gradient criterion. Evaluating 
Equation 5.3 for the largest particles in Table 3.1 and for the breakdown fields given 
there shows that nowhere is the criterion satisfied. The closest approach occurs for 
Case 14 with relative air density .3, for which the gradient is 3.72 x 10 8 V/m 2 . This is 
close enough to the experimental value that it is conceivable that an occasional larger 
particle may exist which could sustain the spark discharge. 
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6.0 


MULTIPLE PARTICLE EFFECTS 


The discussion of the last section was confined to the field distribution 
around single particles. It was found that it is relatively easy to produce breakdown 
fields around a typical thunderstorm particle, but only for the very largest particles is 
the criterion of maximum field gradient satisfied. However, if one allows the possibility 
of interacting particles, the situation becomes much more favorable for spark initiation. 
By interaction of particles is meant the close approach of two or more particles such 
that the field distributions around the particles are altered from what they would be in 
the individual particle case. In principle, by bringing two particles together closely 
enough, it is possible to achieve as low a field gradient as desired, while still 
maintaining breakdown field intensities. This allows the initiation of a spark between 
the particles. The problem is then to determine what happens to that spark when it 
attempts to propagate away from the two particle system. In effect what has been 
created by the spark between the two particles is a larger elongated particle which 
may have properties more conducive to the initiation of a real lightning leader. 

Systems of particles have not been investigated in detail in this study. To do 
so realistically requires the use of statistical descriptions of particle sizes, densities, 
and motions within a thundercloud. For the scenario described above, however, it 
may be possible to apply the analysis of Sections 2 and 3 to the multiple particle 
interactions. It may be possible to derive information on the density of interacting 
systems of particles necessary to produce a reduction in the ambient field needed to 
produce lightning initiation. This, along with measured and calculated data from 
actual thunderstorms, may allow a statistical prediction of the locations for lightning 
initiation in a cloud. 

7.0 SUMMARY 

In this report, the effeci of thunderstorm particles on lightning initiation has 
been studied. The following major points may be summarized. 

(1) A model has been developed to roughly predict breakdown field as a 
function of particle size, shape, and density. 
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(2) The model has been applied to typical thunderstorm particle 
environments. The results show that only the largest particles may be 
expected to have much effect on the threshold initiation field. 

(3) The effect of net charge on the particles is to lower the effective 
breakdown field. It also makes the local corona breakdowns around 
particles one-sided. 

(4) An experimental criterion on electric field gradient has been applied to 
the fields around thunderstorm particles. Once again, only the largest 
particles were seen to satisfy this criterion. 

(5) Although interactions between particles were not investigated in detail, 
it is clear that close approaches between two or more particles more 
easily satisfy the two criteria on field magnitude and gradient for spark 
initiation. 

(6) Statistical approaches are necessary to investigate the likelihood that 
particle systems are responsible for the initiation of lightning 
discharges. 
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